First-principles modeling of ferroelectric capacitors via constrained-D calculations 
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First-principles modeling of ferroelectric capacitors presents several technical challenges, due to 
the coexistence of metallic electrodes, long-range electrostatic forces and short-range interface chem- 
istry. Here we show how these aspects can be efficiently and accurately rationalized by using a 
finite-field density-functional theory formalism in which the fundamental electrical variable is the 
displacement field D. By performing calculations on model Pt/BaTiOs/Pt and Au/BaZrOs/Au 
capacitors we demonstrate how the interface-specific and bulk-specific properties can be identified 
and rigorously separated. Then, we show how the electrical properties of capacitors of arbitrary 
thickness and geometry (symmetric or asymmetric) can be readily reconstructed by using such in- 
formation. Finally, we show how useful observables such as polarization and dielectric, piezoelectric 
and electrostrictive coefficients are easily evaluated as a byproduct of the above procedure. We ap- 
ply this methodology to elucidate the relationship between chemical bonding, Schottky barriers and 
ferroelectric polarization at simple- metal/oxide interfaces. We find that B02-electrode interfaces 
behave analogously to a layer of linear dielectric put in series with a bulk-like perovskite film, while 
a significant non-linear effect occurs at AO-electrode interfaces. 

PACS numbers: 71.15.-m 73.30.-fy 73.61.-r 77.22.-d 77.65.-j 77.80.-e 77.84.-s 
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I. INTRODUCTION 

Capacitors based on ferroelectric perovskites hold 
promise for substantial advances in nanoelectronics, 
with potential applications in non- volatile random-access 
memories and high-permittivity gate dielectricsii Thin- 
ner devices, which are mandatory for optimal efficiency 
and speed, are strongly influenced by the electrical and 
mechanical boundary conditions imposed by the inter- 
facei^ While there has been significant progress in the un- 
derstanding of strain effects,"^ the electrostatics of metal- 
ferroelectric interfaces still remains a challenge, and is 
widely recognized as a central issue in the scaling of fer- 
roelectric devices. 

Interface electrostatics is generally modeled, within 
Landau-Ginzburg theories, by a hypothetical thin layer 
of standard dielectric ( "dead layer" ) interposed between 
an ideal electrode and the active, bulk-like ferroelectric 
film. The dielectric dead layer is arranged in series with 
the film, and therefore the small interfacial capacitance 
associated with it tends to suppress the polarization of 
the film via a depolarizing field.'* It was postulated a long 
time ago^ that, even in the absence of an extrinsic inter- 
facial layer, a small interfacial capacitance can originate 
from the finite penetration length of the electric field in 
a real electrode. The imperfect-screening model and the 
dead-layer model are mathematically equivalent and lead 
to the same consequences, regardless of the microscopic 
nature of the effect 

Owing to the complex structure and chemistry of a re- 
alistic interface, however, it is difficult to infer the mag- 
nitude of this interfacial capacitance based on macro- 
scopic considerations. Moreover, the usual assumption 



that the capacitance (or equivalently, the effective screen- 
ing length) is constant as a function of the ferroelectric 
displacement might not be justified in some cases. For 
example, it was shown very recently by means of first- 
principles calculations that chemical bonding across the 
junction profoundly infiuences the ferroelectric proper- 
ties of the device.^ This is likely to introduce nonlin- 
earities in the electrical response of the interface that 
are neglected within most phenomenological approaches. 
In order to achieve a quantitative model of the elec- 
trode/ferroelectric interface there is therefore the clear 
need for a theory that complements Landau free-energy 
expansions with a microscopically reliable description of 
local chemistry and electrostatics. 

A strategy for modeling the ferroelectric behavior 
of symmetric and asymmetric capacitors that combines 
Landau theory with first-principles calculations was re- 
cently proposed by Gerra et al? Their strategy has the 
advantage of exploiting the power of the ab-initio ap- 
proach to gain quantitative insight into the coefficients 
that describe the behavior of the interface. In particu- 
lar, the interface enters the free energy through two dis- 
tinct quadratic terms, a depolarizing effect which pro- 
vides a uniform electric field and is the main contribu- 
tion, and a short-range chemical bonding effect which 
provides a much smaller correction. These coefficients 
are then input into a standard Landau free-energy expan- 
sion and used to predict the behavior of devices of macro- 
scopic thicknesses, which are not directly tractable from 
first principles. This model was shown to describe the 
SrRuOa/BaTiOs/SrRuOa system quite accurately. The 
results were consistent with the seminal work of Junquera 
and GhosezfSi who demonstrated how the main impact of 
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the electrodes is embodied in the depolarizing field in an 
otherwise bulk-like BaTiOs film. 

Some authors, however, have questioned the general- 
ity of such an assumption, postulating that in some cases 
the electrodes can have a much more profound impact on 
the ferroelectric film. The authors of Ref. [13, for exam- 
ple, claimed that SrRuOa electrodes can destroy the po- 
lar soft mode of ferroelectric KNbOs films, producing a 
head-to-head domain wall a few unit cells from the inter- 
face. Furthermore, in Ref. [ill, Pt electrodes were found 
to induce a "ferrielectric" dipole pattern in the whole vol- 
ume of a BaTiOs film. Such effects, which are nonlocal 
in nature, cannot be described by the simple models of 
RefsH ancQ. To account for (and clarify the nature of) 
these "exceptions" , it would be very desirable to have a 
rigorous methodological framework that treats the elec- 
trical properties of a given capacitor heterostructure fully 
from first-principles, without any a priori assumptions. 

Such a methodological framework was recently devel- 
oped for the case of purely insulating perovskite super- 
lattices. By performing the calculations at a fixed value 
of the electric displacement field D, Wu et aiH were able 
to separate the long-range electrostatic interactions be- 
tween layers from the short-ranged compositional depen- 
dence. Based on this separation, the electrical properties 
of a given layer were shown to depend on the chemical 
nature of a small number of first and second neighbors 
only. This allowed for a first-principles description of 
dielectric, ferroelectric and piezoelectric properties of ar- 
bitrary superlattice sequences in terms of very few pa- 
rameters, appropriately arranged in the form of a cluster 
expansion. 

It is the main scope of this work to extend these ideas 
to the case of ferroelectric films with metallic electrodes. 
Such an extension is now possible, as there are well- 
established methods for treating polarization and electric 
fields in metal/insulator heterostructures, and these can 
readily be combined with recently-developed approaches 
for treating the electric displacement field D as the con- 
trolled electric variable i^^i^^ Using an extensive analy- 
sis of several Pt/BaTiOs/Pt and Au/BaZrOs/Au capac- 
itor heterostructures to illustrate the power of this ap- 
proach, we show that a film-electrode interface behaves 
analogously to an insulator-insulator interface in a fer- 
roelectric superlattice (assuming that there is no Schot- 
tky breakdown), in that the same "locality principle"— 
holds. This means that the film is in a bulk- like state 
except for the two or three oxide monolayers which lie 
closest to the boundary. Moreover, all the complexity 
of interfacial chemical bonding and electrostatics can be 
incorporated in a single energy contribution, which we 
define as the interface electric equation of state. Tak- 
ing advantage of the constrained-!? technique, we fur- 
ther show how to extract in practice (from calculations of 
compositionally symmetric capacitors) such an interface 
equation of state, and represent it in terms of a potential 
drop which is in general a nonlinear function of the elec- 
tric displacement field. Then, we use this information. 



together with the bulk equation of state of the ferroelec- 
tric, to predict, with full first-principles accuracy, the 
electrical properties of capacitors of arbitrary thickness 
and geometry (symmetric or asymmetric). Finally, we 
show how useful observables such as polarization and di- 
electric, piezoelectric and electrostrictive coefficients are 
easily evaluated as a byproduct of the above procedure. 

Our results demonstrate the validity of Z? as a fun- 
damental electrical variable to study ferroelectric capaci- 
tors within an "imperfect screening" regime. (The appro- 
priateness of such an approach was recently questioned, 
although in a slightly different context, in Ref. [13 .) 
From the practical point of view, our detailed study of 
Au/BaZrOa/Au capacitors also yields important insight 
into the similarities and dissimilarities of AO-terminated 
versus B02-terminated perovskite films in contact with 
simple-metal electrodes. On the one hand, the relatively 
high interfacial capacitances we obtain for both interface 
types corroborate the ideas of Ref. 7, where weak inter- 
face bonding was found to be favorable for the overall 
dielectric (or ferroelectric) response of the device. On 
the other hand, at the BaO-Au interface we find signifi- 
cant non-linear effects, which do not fit into a "constant 
interfacial capacitance" model. We correlate these effects 
with the formation and breaking of the interfacial Au-0 
bonds upon polarization reversal. (This same mechanism 
was already found to strongly influence the ferroelectric 
instability in Ref. [3.) 

The manuscript is structured as follows. In Section HIl 
we review the methodological background and present 
the new developments which are specific to this work. In 
Section Uni we discuss the structural and electronic prop- 
erties of ferroelectric Pt/BaTiOa/Pt capacitors, which 
we then use to model their dielectric and piezoelectric 
properties as a function of thickness and applied bias. In 
Section |TV] we focus on the Au/BaZrOa/Au model sys- 
tem. First we compare the electrical properties of the 
Au-BaO and the Au-Zr02 interface structures; then we 
show how to reconstruct the behavior of asymmetric ca- 
pacitor configurations starting from the interfacial and 
bulk equations of state. Finally, in Section |V] and Sec- 
tion IVII we discuss our results in light of the existing 
literature and present our conclusions. 

II. METHODS 
A. Polarization 

1. Bulk insulators 

We shall consider superlattices and capacitor struc- 
tures stacked along z, so that we are interested in po- 
larizations and fields only along this direction. We start 
with the case of a bulk insulator, either a single bulk unit 
cell or a supercell representing an insulating superlattice, 
but with a formulation chosen for convenient later gen- 
eralization to the case of a capacitor structure. 
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We thus consider a periodic insulator described by 
three real-space lattice vectors R^, where for simplic- 
ity of notation we impose that R3 = (0, 0, c) is per- 
pendicular to R1.2 (the latter two lie therefore in the 
xy plane) ; the corresponding reciprocal-space vectors are 
Gi.2,3. We choose a discrete fc-point sampling of the 
form k = j'b|| -|- kx, where the vector b|| = Ga/A^n spans 
a regular one-dimensional mesh of dimension iV|| , and k^ 
belongs to a set of N± special points in the perpendicular 
plane. The electronic ground state is defined by a set of 
one-particle Bloch orbitals, Unk] our goal now is to define 
the polarization along G3. 

To that end, we first seek a localized representation of 
the electronic wavefunctions along the direction G3 for 
each given kj_ . We do this by constructing a set of max- 
imally localized "hermaphrodite" orbitals w„ki (r) that 
are Wannier-like along z while remaining Bloch-like in 
the xy planei^ii^ using a parallel-transport procedure^ii 
The center Znu± of w„kj_ is then defined asiS 

r 

Zn\ij_ = {Wnk^\z\Wnk^) = \Wnk^ ir)f zdr^ , (1) 

and the contribution of k^ to the polarization is 

^'(ki) = ^(-2e^z„k^ +^g„z„), (2) 

n a 

where Za and Qa are the ionic coordinate and bare pseu- 
dopotential charge, respectively, of the atom a; the factor 
of two refers to spin-paired orbitals. The total polariza- 
tion P is then obtained by averaging k_L over its 2D Bril- 
louin zone, while being sure that the branch choice of 
P(k_L) is continuous in k^. 

Our Wannier-based definition of P, Eq. ([2]), lends itself 
naturally to a local decomposition in terms of the dipo- 
lar contribution of individual oxide layers, as proposed 
in RefsH^ an(Jl2l. In particular, given that in typical 
perovskite insulators the centers Znk± cluster themselves 
around the oxide layers they formally "belong" to, one 
can define the layer polarization (LP) of the j-th layer as 

P.(k±)-|(-2^z„k,-|-^Q„Z„). (3) 

In the above equation the sums are restricted to atoms 
a and Wannier centers i that are "in" the layer j, and S 
is the cell surface area; again, the overall pj is calculated 
by performing a 2D Brillouin zone average, pj is well 
defined as long as (i) the oxide layers are charge-neutral, 
and (ii) the assignment of a specific atom or wavefunc- 
tion to a given layer is clear-cut and unambiguous. Both 
conditions are satisfied in typical II-IV perovskite ferro- 
electrics such as PbTiOs and BaTiOa. 

2. Capacitor superlattices 

Ideally one would like to study a capacitor in the form 
of a number of layers of insulator sandwiched between 



semi-infinite metallic contacts. However, we adopt here 
the standard approach of constructing supercells consist- 
ing of alternating insulating and metallic regions stacked 
along z, just as is normally done when studying ferro- 
electric superlattices. We adopt the same notations and 
conventions as in the previous subsection, with c being 
the superlattice repeat distance along z. We set A^|| = 1 
(and henceforth write k_L = k); this is by no means a lim- 
itation, since we are only interested in capacitors that are 
thick enough so that tunneling is insignificant, in which 
case the one-particle bands will have negligible disper- 
sion along the z direction. We further require a rectify- 
ing (rather than ohmic) contact at the oxide/electrode 
interface. This means that both the valence-band max- 
imum (VBM) and conduction-band minimum (CBM) of 
the film are located far enough in energy from the Fermi 
level that they are not appreciably populated/depleted 
by the tails of the smearing function (e.g., Fermi-Dirac, 
Gaussian, etc.). 

Because the capacitor superlattice is metallic, one 
might wonder whether it is possible to define a polar- 
ization P. However, the superlattice is only metallic in 
the X and y directions, whereas we are interested only in 
computing P along z, and only in applying fields along 
z. The methodology for computing P in such cases was 
developed in Ref. [20. The electronic states are classified 
into three energy windows: 

• The completely empty states (upper window) are 
discarded from the computation, since they do not 
contribute to P or to other ground-state properties. 

• The partially occupied states (middle window) ly- 
ing in the range W = [Ep — 6,Ei + S] around 
the Fermi level E{ are considered as conduction 
states. Since these states fall in the energy gap of 
the dielectric film, they are confined to the metallic 
slab, and the dipole moment of their overall charge 
distribution is thus well defined. To make sure 
that this conduction charge distribution decays fast 
enough in the insulating film, it is useful to define 
its planar average 

Pcond{z)^^ ^ Wk/«k / rfa;dy |7/;„k(r)P , (4) 

where e„k and /„k are the eigenvalue and occu- 
pancy of the state, Wk is the fc-point weight, and S 
is the cell cross-sectional area. 

• The lower states, which are all fully occupied, are 
transformed to yield a set of hybrid Wannier func- 
tions Wkn that are maximally localized along z 
while remaining extended (and labeled by k = kx) 
in the x and y directions. The contribution of 
each Wannier function to P is then defined through 
the center of the corresponding charge distribution 
Pkn{x) = |u;kn(a;)P- 

The center of charge of Pcondiz) (middle window) is com- 
puted by integrating against a linear sawtooth function 
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whose discontinuity is placed in the middle of the in- 
sulating region. Similarly, the center of each Wannier 
charge (lower window) is computed from its pkn using a 
sawtooth function whose discontinuity is chosen far away 
from its center. For the systems considered in this work, 
the pkn SLie typically very well localized, and the main 
source of error comes from the slower decay length of pc 
in the insulator. This means that in very thin capacitors 
(few oxide layers) the polarization becomes ill-defined; 
rather than a defect of the algorithm, this is a signature 
that the system becomes metallic, and the polarization 
cannot be defined. 

Note that there is an inherent arbitrariness in the sep- 
aration of the total charge density into lower and middle 
windows (i.e., in the choice of the parameter S above). 
This arbitrariness indeed affects both pcond and those 
Pkn which lie closest to the electrode; the total value of 
P, however, is not affected by this choice, and is there- 
fore well defined. Far enough from the electrode, the pkn 
themselves are unaffected by this choice, and can there- 
fore be used to construct meaningful layer polarizations, 
analogously to the case of an insulating superlattice. This 
point will be demonstrated in practice in the applications 
sections. 

We are generally concerned with capacitor structures 
in which a finite bias is applied across the capacitor. In 
our approach, this is treated by applying a finite macro- 
scopic electric field £ along the z direction of the super- 
lattice, and identifying £c as the bias applied between 
successive metallic segments. The formulation above ap- 
plies equally well to this case, where it is understood that 
the electric field couples to pcond (middle window) and to 
all the Wannier charges (lower window). Note that the 
presence of a finite macroscopic field implies that there is 
effectively an infinite number of regularly spaced Fermi 
levels, one for each repeated image of the metallic slab 
along the field direction. The "transition" between two 
adjacent Fermi levels takes place deep in the insulating 
slab, where the system is locally insulating and a shift in 
E{ within the gap does not produce any physical conse- 
quence. 



B. Constrained-D method 

1. General theory 

We summarize here the details of the constrained 
displacement-field method that are most relevant for this 
work (see Ref. J^for the full derivation). For consistency 
with the previous sections, we restrict our analysis to the 
case of a monoclinic system, with the polarization axis, z, 
parallel to the heterostructure stacking direction and per- 
pendicular to the xy plane; we shall further assume that 
Ri^2 are fixed, and only c (together with the ionic and 
electronic coordinates, {v}) is allowed to vary. Within 
these assumptions, the constrained-D method^^ reduces 
to a simpler formulation, where only the z components of 



the macroscopic fields D, P and £ are explicitly treated. 
Thus, we define the internal energy functional 

U{D, {v}, c) = Eks{{v}, c) + ^[D^A^P{{v}, c)\\ (5) 

which depends directly on the external parameter £), 
and indirectly on the internal ({w}) and strain (c) vari- 
ables through the Kohn-Sham total energy Eks and the 
macroscopic polarization P; S ~ |Ri x R2I is the con- 
stant cell cross-section. We then proceed to minimize the 
functional with respect to v and c at fixed D: 



U{D) = min [/(£), {w},c). 



(6) 



which yields the equilibrium state of the system as a func- 
tion of the electric displacement D. 

D can also be expressed in terms of the reduced vari- 
able d = SD/Att, which has the dimension of a charge 
and can be interpreted as d = —Qhee, where Qfioc is the 
free charge per surface unit cell stored at a hypothetical 
electrode located at z = +00^ Since the surface areas 
of the parallel plate capacitors considered in this study 
are not allowed to vary, constraining D or d is completely 
equivalent. However, for reasons of convenience, we shall 
adopt d as our electrical variable in the remainder of this 
work. 

This method is equally valid for bulk insulators, in- 
sulating superlattices, and capacitor superlattices, once 
the polarization is defined as explained in Sec. Ill Al For 
the capacitor case, our adoption of the definitions of 
Sec. Ill A 21 implies that the metallic electrode layer is 
treated as an infinitely polarizable dielectric, and the free 
charges on its surfaces are reinterpreted as polarization 
charges coming from the metal. While such a choice may 
seem unnatural from the point of view of textbook elec- 
trostatics, it is in fact the most natural one in the con- 
text of first-principles electronic-structure calculations, 
where it is not easy to draw a distinction between free 
and bound charges. For example, the metal-insulator in- 
terface is typically rather diffuse, with the conduction 
states of the metal mixing strongly with the states of the 
insulator across several interatomic spacings, so that a 
spatial distinction is not meaningful, and we have seen 
in Sec. Ill A~2l that a distinction based on energy windows 
is also arbitrary to some degree. 

The reduced electric field e = £c, which is minus the 
potential step across the supercell, e = —V, is related to 
the internal energy by 



e{d) 



dU{d) 
dd ■ 



(7) 



This corresponds to the fundamental relationship 



U{D2) - U{D,) = ^ £{D)dD. (8) 



Di 



of classical electrostatics, but expressed in differential 
form using the reduced variables appropriate to the 
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variable-cell case. The connection to classical electrostat- 
ics can be made even more apparent by recalling the rela- 
tionship between the reduced variables and free charges 
and potentials 



D 



U[D2) - U{Di] 



d2 



e{d)dd 



V{Q)dQ. (9) 



Having established the functional relationships be- 
tween the active degrees of freedom (both electrical and 
structural), it is relatively easy now to extract from a 
calculation all functional properties of a material or de- 
vice that involve a coupling between them. For example, 
the proper piezoelectric strain constant can be readily 
obtained as 



-^33 



dc dc ^ de^~^ 
de dd \dd/ 



(10) 



Note that in the above equation de/dd has the dimension 
of an inverse capacitance and is related to the free-stress 
dielectric constant of the crystal; with the notation of 
Ref. \2fA we have 



^33 



47rc /de\-i 
~\dd) ' 



(11) 




FIG. 1: (Color online) Sktech showing conservation of longi- 
tudinal component of displacement field D, but not electric 
field £ or polarization P, in an insulating superlattice com- 
posed of three dielectric constituents A, B and C. 



C. Locality principle and spatial decomposition 

According to classical electrostatics, in the absence of 
free charge the normal component of the electric displace- 
ment field is preserved at a planar interface between two 
insulators. 



(D2-Di)-n = 0. 



(12) 



2. Practical procedure 

We typically span the range of relevant polarization 
states by repeating the structural and electronic relax- 
ations for a number (five to ten) of equally spaced d 
values, d = di, ^2, c?„. For each di, we tabulate the 
energy C/^, (cell-averaged) electric field £i, and equilib- 
rium out-of-plane lattice constant c^; we use the latter 
two to compute the reduced field — £iCi. 

In order to obtain Si with sufficient accuracy it is im- 
portant to perform well-converged structural relaxations, 
by imposing sufficiently stringent thresholds on residual 
forces and stresses. A useful indicator of the overall 
numerical quality of the calculations is the fundamen- 
tal relationship Eq. ([7]), which in principle should be 
exactly satisfied. To check the validity of Eq. ([7]) we 
first perform a polynomial fit to the calculated {ei,di) 
points. This yields a continuous function, e{d) that we 
integrate analytically to obtain U{d) modulo a constant; 
we choose this constant as the one that best matches 
the first-principles internal energies {Ui, di). Any residual 
discrepancy between Ui and U{di) points to a numerical 
issue that must be addressed before proceeding further in 
the analysis. Usually the most important source of error 
concerns the relaxation of the cell volume and shape; we 
shall discuss this issue further in Section UTeI 

For convenience, we also perform a polynomial fit to 
the {ci,di) points, which yields a continuous curve c{d) 
that is relevant for the piezoelectric response of the crys- 
tal, as stated in the previous subsection. 



This means that, for an insulating superlattice in elec- 
trostatic equilibrium, D is the same in all participating 
layers, unlike the electric field £ and the polarization P, 
whose local values generally vary from layer to layer (see 
Fig. [1]). Therefore, using D (or d) as the fundamental 
electrical variable is extremely practical for modeling the 
behavior of ferroelectric capacitors, because it makes it 
possible to decompose the equation of state of a layered 
structure into the sum of the individual building blocks. 
For example, we can write the internal energy as 

Uid)^Y.Wl (13) 

i 

where Ui refers to the internal energy of an appropri- 
ately defined sub-unit. For a capacitor with metallic elec- 
trodes, it is natural to decompose the internal energy as 

Uid) ^Ui^ + Ui^{d) + 7V[/b(d) + C/R(d) + Un (14) 

where N is the number of bulk cells comprising the in- 
sulating film and Uh is its bulk internal energy per cell, 
Uj^{d) and U^iid) are the left (L) and right (R) interface 
internal energies, and C/l and C/r are the internal energies 
of the left and right metallic electrodes. (In our capacitor 
supercells, I7l and J/r, are combined into A^mctait^mctai, 
where iVmetai is the number of cells of bulk metal and 
C^mctai is its internal energy per cell, which is independent 
of d as appropriate for a metal.) Taking the derivative of 
Eq. pi]) according to Eq. ^ yields 

e{d) = si^id) + iVeb(d) + £R(d), (15) 
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where Ebuik is the potential drop across a unit ceU of the 
bulk insulator at a given value of d, and £l.r — dUi^^^/dd 
contains the interface-specific information. 

The potentials and the energies contain the same infor- 
mation, apart from a constant of integration, and one can 
choose to work with one or the other as a matter of prac- 
tical convenience. Indeed, when analyzing the electrical 
properties of a capacitor, one is generally interested in 
energy differences between two different electrical states, 
rather than in the total energy of the device. There- 
fore, the constant of integration that gets lost in going 
from Eq. to (|15p is not important for the scope of 
our discussion. Thus, we shall assume henceforth that 
t/(0) = 0, which also implies that the constant energies 
C/l,r in Eq. (I14|) have been set to zero. 
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FIG. 2: Schematic model of the decomposition of the poten- 
tial into bulk and interface contributions. An electron trav- 
eling from the left (L) electrode to the right (R) electrode 
experiences potential variations of —<j>h, A''£buik and (j>R\ the 
total variation is £ = -Ef (R) — -Ef (L). 



D. Decomposition of the interface contribution 

1. Partial decomposition 

Since U{d) and e{d) in Eqs. (|14m5p can be obtained 
from supercell calculations, while Uh{d) and eh{d) can be 
obtained from bulk insulator calculations, it is straight- 
forward to extract the quantities 

C/i„t(rf)-C/L(d) + C/R(d) (16) 

and 

eint(d) =eL(d) + £R(d) (17) 

representing the total impact of both electrodes on the 
electrical equation of state of the capacitor. Explicitly, 
we take 

eint{d) = eN{d)-Neb{d). (18) 

Often, this is all that is needed, e.g., for modeling the 
polarization and dielectric response of a given device as 
a function of the oxide film thickness (we shall demon- 
strate this in our first application to Pt/BaTiOs/Pt ca- 
pacitors). The number N of cells of insulating material 
should be kept small enough to avoid an undue computa- 
tional burden, while remaining large enough to decouple 
the two electrode interfaces, so that the center of the 
oxide slab should behave like the bulk material within 
the same mechanical (in-plane strain) and electrical (d) 
boundary conditions. 

2. Full decomposition 

There are, however, situations in which it may be valu- 
able to obtain the individual terms in Eq. ([TT]), i.e., to 
define the individual interfacial potential steps £l and er 
which occur at the left and right electrode interfaces re- 
spectively. However, instead of using quantities defined 
as offsets of the average electrostatic potential across the 



interface, we find it more physical to use variables (^l and 
(/>R that are the offsets of the metal Fermi levels Ef{L) 
and Ef{R) relative to the VBM just inside the insula- 
tor, as illustrated in Fig. (5] With this choice, and (j)K 
are just the p-type Schottky barrier heights (SBH) at the 
metal/insulator interface. (It would be equally viable to 
adopt the CBM as the reference, corresponding to n-type 
Schottky barriers, but we do not do so here.) As long as 
both electrodes are made from the same material;^ the 
total potential step is just e = Ei{R) — E[{L). This dif- 
ference can be decomposed by following the hypothetical 
path in Fig. [5] of an electron traveling from the left elec- 
trode through the insulator and into the right electrode, 
and we obtain 

e = ^<j)L+Neb + (t)n. (19) 

In the remainder of this section, we make some of the 
above definitions more precise, and discuss how in prac- 
tice to extract accurate values of the SBH at a polar- 
ized metal/insulator interface. The main issue here is 
that, whenever Sh is non-zero, the SBH is somewhat ill- 
defined because the VBM does not have a well-defined 
asymptotic value deep in the oxide. (Instead, it varies 
linearly with depth, with a slope corresponding to the 
internal electric field £h)- In the next few paragraphs, 
we propose a procedure that provides a reasonable yet 
sharp definition of and even when 5b 7^ 0. 

For the moment we assume an interface configuration 
with the semi-infinite electrode at right and the film at 
left. The situation is sketched in Fig. [31 which also il- 
lustrates the following discussion. We first compute the 
planar average of the local electrostatic potential, Vf/(r), 
further convoluted with a Gaussian filter of width a to 
suppress the short-range oscillations, 

Vid, z) = [ Vnid, r')e-(^-^')'/"'d3r' . (20) 

Next, we identify two z coordinates on either side of the 
interface, zp in the film and zm in the metal. Both zp 
and Zm must be located far enough from the interface 
that the short-range structural distortions related to in- 
terface bonding have already relaxed back to the regular 
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FIG. 3: (Color online) Illustration of proposed procedure 
for evaluating the p-type Schottky-barrier height at a per- 
ovskite/metal interface when a macroscopic field is present in 
the insulating layer. Light blue circles correspond to B-site 
cations, red circles to oxygens, and gold circles to the metal 
electrode atoms. 



bulk-like spacings of the respective material (oxide film 
or metal electrode). As a further requirement, we im- 
pose that the interface-related perturbations in the local 
electrostatic potential V{z) (schematically indicated in 
the figure by the strong oscillations near the metal/film 
boundary) are also negligible in the neighborhood of both 
Zp and Zm- This implies that V{d, z) is a linear function 
near zp with a finite slope given by the bulk internal field 
fbuik('i), and V{d,z) is a constant near zm- We gener- 
ally find that three perovskite unit cells on the film side 
and five monolayers on the electrode side are sufficient 
for both requirements (on the structure and V{z)) to be 
satisfied accurately. We therefore set zp a.s the z coordi- 
nate of the 4-th B-site cation in the film (numbering as 
1 the B cation which lies adjacent to the interface), and 
Zm as the z coordinate of the 6-th metallic layer in the 
electrode. 

Now, using these two reference points in the lattice, 
we extract V{d,ZM) and V{d,zp), which are indicated 
in the figure as black circles. On the film side, we use 
V{d, Zp) to estimate the VBM at the interface, 



VBM 



(d) = V{d, zp) + AVp{d) + 3£-buik(d), (21) 



where AVp{d) is the relative position of the VBM to the 
average electrostatic potential in the bulk (techniques for 
calculating this are detailed in the following subsection). 
On the metal side we compute the Fermi energy as 



Ef{d)^V{d,ZM) + AVM, 



(22) 



where AVm (independent of d) is again a bulk property, 
i.e., the relative position of the Fermi level to the average 
electrostatic potential of the metal. Finally, we define the 
p-type Schottky barrier as 



It is easy to verify that this definition reduces to the 
standard technique for calculating Schottky barriers at 
metal/semiconductor interfaces^! whenever the macro- 
scopic field in the oxide vanishes. Note that the above 
construction provides, as a byproduct, structural pa- 
rameters that are relevant for accessing the piezoelectric 
properties of the device. In particular, starting from the 
same zp and zm, we define an interfacial expansion 



c{d) = \zM - zp\ - iVcbuik(rf) - 5(5j\/, 



(24) 



^id)=Et{d)-EyBM{d). 



(23) 



where Sm is the bulk interlayer distance of the metal (see 

Fig. ED. 

Note that, while there is an intrinsic arbitrariness in 
the definition of (f>{d) and c{d) (several choices are pos- 
sible for Zp), the arbitrariness always cancels out in the 
final equation of state of the entire capacitor because 
of the way these functions are always summed up in 
pairs (a capacitor always has two electrodes). We also 
note that these functions, by construction, transform 
properly under spatial inversion, so that for a capacitor 
having a centrosymmetric reference structure, we have 
(I^Rid) 4>L{-d) and CR{d) = CL{~d). 



3. AVf and AVm 

AVm can be calculated with high precision for the 
bulk metal by extracting the Fermi level and the aver- 
age electrostatic potential from the structural and elec- 
tronic ground state. To define AVp{d) we start from 
a constrained-D calculation of the bulk oxide. Since a 
macroscopic electric field is generally present, the values 
of both the VBM and the average electrostatic potential 
are not directly obvious from the eigenvalue spectrum 
(strictly speaking, the energy eigenvalues themselves are 
ill-defined). The effect of an electric field is to induce 
a linear ramp in the electrostatic potential, and a cor- 
responding linear "tilting" of the energy bands. For a 
given value of the macroscopic electric displacement, the 
VBM and the average electrostatic potential will there- 
fore have the same linear z-dependence, and the differ- 
ence AVf = Vyb.m{z) — Vh(z) will be independent of z. 
In practice we compute AVp by first relaxing the struc- 
tural and electronic degrees of freedom in the finite field, 
using the usual convention that the electrons feel a peri- 
odic electrostatic potential having zero unit-cell average, 
plus a coupling to the field through the Berry-phase po- 
larization. AVp is then obtained by diagonalizing the 
zero-field Hamiltonian operator in the subspace spanned 
by the wavefunctions, which form the "ground state" of 
the finite-field calculation, and finding its maximum over 
the wavevectors in the Brillouin zone. 

Note that this procedure is to some extent arbitrary, 
and it is certainly possible to adopt alternative strategies. 
Whatever choice is made, the only important require- 
ment is to have a well-defined reference energy in the 
insulating lattice as a function of D] the arbitrariness in 
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the specifics of tfiis choice cancel out anyway when we 
consider a complete capacitor heterostructure. 



E. Computational parameters 

Our calculations are performed within the local- 
density approximation of density-functional theory and 
the projector-augmented-wave method^^ as implemented 
in an "in-house" code. We used a planewave basis cutoff 
energy of 40 Ry in Section |TTT] and of 80 Ry in Section 
IIVI the higher value in the latter case is intended to min- 
imize the Pulay error in the stress and the numerical 
noise in the energies which are due to the discrete na- 
ture of the basis seti^ In all cases we fix the in-plane 
lattice parameter to a constant value and we enforce a 
tetragonal PAmm symmetry constraint; the out-of-plane 
lattice parameter, as well as the internal coordinates, are 
allowed to relax fully. The Brillouin zone integrations 
of the capacitor heterostructures are performed with a 
6x6x1 mesh, where kz — and the grid is shifted 
in-plane according to the Monkhorst-Pacl^ prescription 
for two-dimensional sampling; the Gaussian smearing en- 
ergy is set to 0.15 eV. In the bulk calculations we use a 
6x6x6 Monkhorst-Pack mesh, which is sufficient to 
converge both the structural and the dielectric response 
of the crystal to an accuracy comparable to that of the 
capacitor calculations. To relax the structure (both in- 
ternal coordinates and the out-of-plane strain) at each d 
value we use a steepest-descent approach, optimally pre- 
conditioned by inverting the force-constant matrix and 
the elastic constant calculated in the centrosymmetric 
d=0 configuration. Generally, five to ten iterations were 
sufficient to relax the geometries to a stringent conver- 
gence threshold for both forces (10~^eV/A) and stresses 
(10 MPa). (To ensure excellent accuracy of the calculated 
energies and potentials, we further enforce a threshold of 
1% convergence in the internal electric field.) 

Correcting for the Pulay error in the stress is crucial to 
accurately model the strain-polarization coupling effects 
discussed in this work. We use a technique similar in 
spirit to the prescription of Ref. [26l . In particular, we 
define the corrected stress (T^j as 



where (t°- is the calculated stress tensor (analytical 
derivative at fixed number of plane waves), is the cell 
volume, and C is a constant (dependent on the cell stoi- 
chiometry and plane- wave cut off) . To evaluate C several 
techniques are possible. A possible strategy is to fit the 
dependence of the total energy on the plane- wave cut off, 
as discussed in Ref. [2^. In our calculations, we infer C by 
imposing aij = in Eq. (j25[) for a particular configura- 
tion of a given system that has been structurally relaxed 
using a different technique (e.g., a Murnaghan fit to the 
energy/volume curve). 



III. RESULTS: Pt/BaTiOg/Pt CAPACITORS 

A. Motivation 

Our goals in this section are threefold. First, we shall 
introduce our methods for computing the macroscopic 
polarization in short-circuited ferroelectric capacitors, 
separating the different contributions that we discussed 
in Section [il A 21 Second, we shall analyze the structural 
and electrical properties as a function of the thickness 
of the short-circuited film, identifying those aspects that 
are common to ferroelectric single- crystal BaTiOa, and 
those that depart from the bulk behavior. Third, we 
shall demonstrate with a quantitative model that even 
these thickness-dependent perturbations can be under- 
stood in terms of the bulk properties of BaTiOa, once 
the interface contribution is properly taken into account. 
This analysis is primarily aimed at verifying in practice 
our "locality principle," which allows us to separate the 
equation of state of a capacitor into two interface contri- 
butions and a bulk- like term. We shall show that, in the 
case of Pt/BaTiOa/Pt, this separation works with excel- 
lent accuracy down to a thickness of only two BaTiOs 
unit cells. 

The choice of Pt and BaTiOa, and more specifically of 
the BaO-terminated interface, is motivated by the recent 
prediction^, of a chemical bonding mechanism that en- 
hances the ferroelectricity of the film beyond the bulk 
BaTiOa value. Because of this effect, it was found 
that Pt/BaTiOs/Pt capacitors remain ferroelectric down 
to a single unit cell of BaTiOa, i.e., there is no criti- 
cal thickness below which the polar instability is sup- 
pressed. Given the practical interest in overcoming the 
usually deleterious size effects in ferroelectric devices, the 
Pt/BaTiOs/Pt system is therefore an appealing test case 
for the present study. We warn the reader, however, that 
the above-mentioned features of the Pt/BaTiOa/Pt sys- 
tem are to some extent anomalous, i.e., they depart from 
the usual understanding of depolarizing effects in thin- 
film ferroelectrics. For this reason, the results presented 
in this section should not be understood as an example 
of the most typical ferroelectric capacitor. Instead, this 
application to Pt/BaTiOs/Pt illustrates how the general 
strategy developed in this work (free from apriori as- 
sumptions) is particularly effective at capturing the pe- 
culiar physics of a highly non-standard case. The Ti02- 
terminated interface of BaTiOa with Pt would perhaps 
have provided a more "regular" example, which in princi- 
ple could have allowed us to trace a closer link with earlier 
first-principles and phenomenological results. However, 
this system is inappropriate because it suffers from the 
band- alignment issues mentioned in Ref.fl In particular, 
we find that the Ti02-terminated BaTiOs/Pt interface 
has charge-spillage problems already when the capacitor 
is in the paraelectric reference structure, thwarting at- 
tempts at defining a polarization or even introducing an 
external bias potential. In Section [IVI we consider a dif- 
ferent ferroelectric/electrode combination (BaZrOs/Au) 
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FIG. 4: (Color online) Layer rumplings, defined as cation 
displacements relative to oxygens, for oxide layers in relaxed 
short-circuited Pt/BaTiOs/Pt capacitors containing 7 (cir- 
cles), 5 (squares), 3 (diamonds), 2 (left triangles), or 1 (up 
triangles) perovskite unit cells. Odd and even layer numbers 
refer to BaO and Ti02 layers respectively, with Layer 1 be- 
ing the BaO layer that is chemically bonded to the Pt. Bulk 
values are shown for comparison as filled symbols connected 
by dashed lines. The polarization is along +z. 

whose B02-type interface is free from band-alignment 
problems; in that case we are able to compare two differ- 
ent interface types and discuss the differences between a 
"standard" and a "non-standard" case. 



B. Structural and dielectric properties at zero bias 

We start by analyzing the polar ground state at zero 
bias of 7, 5, 3, 2 and 1-unit cell thick BaTiOa films with 
compositionally symmetric (the overall spatial symme- 
try is broken upon FE off-centering) BaO terminations 
and Pt electrodes. (A^-unit cell capacitors are actually 
"N + 1/2" perovskite cells thick; e.g., 7V=1 means Pt- 
BaO-Ti02-BaO-Pt.) The Pt electrodes are modeled by 
a 9-layer Pt slab in periodic boundary conditions. We 
fix the in-plane lattice constant to ao = 7. 276 a. u., the 
theoretical equilibrium value for cubic SrTiOs, and we 
allow the out-of-plane lattice parameter of the tetrago- 
nal supercell, as well as the internal coordinates, to relax 
fully. We shall present our results starting from a com- 
parative analysis of the relaxed atomic positions in our 
short-circuited Pt/BaTiOs/Pt capacitors; then we shall 
gradually introduce the ingredients that enter the defini- 
tion of the polarization and its coupling to an external 
field. 



1. Structural properties 

We plot in Fig. 2] our calculated results for the layer 
rumplings of the relaxed capacitors at zero bias, defined 



as the cation displacement relative to the oxygens in the 
same oxide layer; in the same figure we report the calcu- 
lated rumpling values for bulk BaTiOa as a comparison. 
The most striking feature at all thicknesses is the strong 
bucking of layer 1, which is the BaO layer directly in 
contact with the Pt surface on the negatively polarized 
end of the films. (The BaO buckling is also enhanced at 
the positively polarized end, but the effect there is signif- 
icantly smaller.) Quite interestingly, the rumplings of all 
the films are systematically larger than the bulk values; 
furthermore, the enhancement in the structural distor- 
tions becomes more important in thinner films. This is 
unexpected, as the depolarizing effect is known to sup- 
press polarization and symmetry-breaking distortions in 
the ultrathin limit. The mechanism leading to such an 
enhancement is related to the interfacial chemical bond- 
ing effect discussed in Ref. 0; we shall clarify this point 
in the following. 

In Fig. Owe plot the interlayer distances between ions 
belonging to the same sublattice, focusing here on the 3, 
5 and 7-cell thick capacitors only. (Atoms are grouped in 
different sublattices according to their chemical identity. 
Pt and O atoms are further split into three and two sub- 
lattices, respectively, as shown schematically in the right 
panel of Fig.O) The atomic layers closest to the interface 
undergo strong distortions, both on the electrode and on 
the insulator side. The largest effects are located again 
on the negatively polarized end of the film. Here, the 
surface Pt atoms and the BaO ions strongly buckle, with 
the overall effect of reducing the Pt-0 distance (which 
ranges from 2.01 A to 2.04 A in the capacitors consid- 
ered) and increasing the Ba-Pt distance ('^ 2.9 A). These 
features are consistent with the oxygen binding chemi- 
cally to the Pt surface, while the Ba atom repels the Pt 
atom that lies directly underneath. Such a picture was 
proposed, from an analysis of the centrosymmetric ref- 
erence structure, in Ref. 0; here we can see its impact 
on the properties of the fully polarized state of the film. 
At the positively polarized end of the film, the structural 
distortions of the Pt surface are relatively minor, and the 
oxide film does not appear to be chemically interacting 
with the electrode; the Pt-0 distances in all capacitors 
are larger than 3.3 A, and the metal-oxide bonding ap- 
pears to be of purely electrostatic nature. Two mono- 
layers away from the interface, the interlayer distances 
of the BTO film converge to a uniform value, which can 
be understood as the relaxed strain state of the film in 
the capacitor heterostructure. In all cases this value is 
larger than in the equilibrium value of the strained bulk, 
which is indicated in the same figure as a dashed hori- 
zontal line. (The bulk out-of-plane strain was calculated 
by imposing the same in-plane strain as in the capacitor 
calculations; therefore, the effect shown in Fig.[5]is not of 
mechanical origin.) Remarkably, the tetragonality of the 
film increases for thinner capacitors. Since ferroelectrics 
have a strong coupling between polarization and strain, 
this provides additional evidence to the enhancement of 
polarity we already pointed out earlier while discussing 
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FIG. 5: (Color online) Relaxed interlayer distances of short-circuited Pt/BaTiOs/Pt capacitors with oxide thickness of 7, 5 
and 3 unit cells and 9 Pt layers. Black (left) and colored (right) symbols/lines correspond to Pt and oxide layers respectively. 
Vertical axis indicates distances between neighboring cations belonging to the same sublattice (see right panel); horizontal axis 
is the mid-point coordinate. Dashed line indicates the calculated equilibrium out-of-plane lattice parameter of bulk BaTiOa 
strained to the STO in-plane lattice constant. The code for the symbols and colors is schematically explained in the right 
panel, where the relaxed structure of the bottom interface (located at a; ~ 10 A) is shown. Each color corresponds to a different 
sublattice: Ptl (light grey), Pt2 (white), Pt3 (dark grey), Ba (red), Ti (light green), 01 (blue), and 02 (dark green). 



the layer rumplings. 

Note that such a strong coupling makes the results 
very sensitive to the accuracy in the relaxation of the 
out-of-plane lattice constant. To this end, it is crucial to 
properly take into account the effect of the Pulay stress, 
as explained in Section III El In order to check that this 
procedure was effective, we monitored in all capacitors 
the interlayer distance in the center of the Pt slab, i.e., 
the black circle at a; = in the three panels of Fig. [5l The 
maximum deviation in this value was less than 10"'^ A, 
confirming that our structural relaxations are very well 
converged. 

In the next subsection we shall investigate the electri- 
cal properties of the Pt/BTO/Pt capacitors, and demon- 
strate the ferroelectric nature of the enhanced structural 
distortions discussed above. 



2. Polarization and electrical properties 

Our goal now is to evaluate the macroscopic polariza- 
tion of the capacitor heterostructures discussed in the 
previous section. Since all the capacitors are relaxed 
within standard short-circuit electrical boundary condi- 
tions, the macroscopic electric field is zero and the polar- 
ization is equal to the electric displacement field D. 

First we assess the level of accuracy we can expect 
for the value of the polarization computed according to 
the technique discussed in the methods section. To that 
end, we analyze the planar-averaged conduction charge 



Pcond(a^), as defined in Eq. (g]) by setting the width of the 
middle energy window S ~ 0.5 eV. Given the Gaussian 
smearing of cr = 0.15 eV, [Ei — 6,Ei + S] encompasses all 
partially occupied states within an occupancy threshold 
of 10~^. In other words, at the bottom of the window the 
smearing function f{E{ — S) is equal to 1 — 10^^, while 
at the top f{Ei + 5) = 10"^. Recall that all states lying 
higher in energy are discarded; all states lying below this 
window are treated as fully occupied and transformed 
into Wannier functions, as we shall discuss shortly. 

To define the dipolc moment of Pcond(a^) it is essential 
that this function be small in the middle of the insulat- 
ing region. Whenever the film is not thick enough for 
Pcondla^) to decay to zero, this introduces an error in the 
definition of the macroscopic P which can be roughly 
estimated by 

AP ^ pcondiL/2)Sx. (26) 

Here L /2 corresponds to the center of the insulating film, 
and Sx is a length that takes into account the arbitrari- 
ness in the positioning of the discontinuity of the saw- 
tooth function used to define the dipole moment of Pcond- 
We shall assume Sx to be half an oxide layer, or 2 bohr 
units. 

To give an idea of how the conduction charge decays in 
the oxide films, we plot in Fig.[S]the calculated Pcond(2^) in 
the five structures considered. While the thinnest capac- 
itor (1-cell thick, dark green curve) is clearly metallic, in 
all other structures the conduction charge decays almost 
to zero in the central oxide layer. The estimated accuracy 
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FIG. 6: (Color online) Conduction charge density pcond(a;) 
(Eq.gl) for five Pt/BaTiOg/Pt capacitors of Fig. g] 



using the above formula is of the order of 0.5/iC/cni^ 
for the 2-cell thick capacitor, decreases by roughly one 
order of magnitude in the 3-cell thick capacitor, and is 
essentially zero in the 5- and 7-cell structures. In the re- 
mainder of this section, therefore, we shall drop the 1-cell 
structure from our discussion and focus on the remain- 
ing four structures, where the value of the macroscopic 
P can be accurately defined. 

The fully occupied states are transformed, separately 
for each fc-point, by means of the parallel transport algo- 
rithmJi This yields a set of orthonormal orbitals which 
sum up to the same charge density and are maximally 
localized along the polarization direction (their Bloch- 
like character is preserved in plane). Note that the ac- 
tual number of states differs for each fc-point. Therefore, 
when performing the Brillouin zone averages of the po- 
larization, particular care must be taken in order not to 
introduce by mistake a fraction of the quantum of polar- 
ization into the final value. In order to ensure that this 
issue is properly taken care of, it is useful to analyze how 
the Wannier centers distribute in space for each fc-point. 
Upon visual inspection, we find that the Wannier centers 
are characterized by a significant degree of disorder in the 
metallic region, as might be expected by recalling that 
the band structure of the metallic slab is not constituted 
by full energy bands, and in our algorithm it is abruptly 
"cut" at E{ — S. Conversely, in the insulating region, we 
find that the Wannier centers cluster nicely around the 
oxide layers, analogously to what happens in purely in- 
sulating superlattices"'^^; to demonstrate such a behavior 
we plot in Fig. [7] the calculated Wannier centers for the 
2-cell capacitor. The total number of orbitals shown in 
Fig. [7| for each fc-point matches exactly the "nominal" 
number of valence orbitals of the oxide ions. This means 
that the oxide film can be identified as a charge-neutral 
and spatially confined subsystem, whose dipole moment 
can be computed with high accuracy, and potential is- 
sues with the quantum of polarization are therefore com- 
pletely avoided. (Note the increased fc-space dispersion 
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FIG. 7: Wannier centers in the 2-unit-cell BaTiOg film as a 
function of in-plane fc-point in the irreducible 2D Brillouin 
zone, labeled as {i,j) according to ky — {i + 1/2, j + l/2)/6. 
Second and fourth groups of centers correspond to Ti02 lay- 
ers; others are BaO, of which the first and fifth are in contact 
with Pt. P points up. 



in the Wannier centers associated with the bottom BaO 
layer; this is due to the strong perturbation induced by 
chemical bonding with Pt.) 

By combining the ingredients discussed in the above 
paragraphs, we now compute the electric displacement 
D of the films, which is plotted as a function of film 
thickness in Fig. [51 In the thickest 7-layer film D = 
46.1/iC/cm^, which is 17% larger than the spontaneous 
polarization of bulk BTO within the same mechanical 
boundary conditions (shown as a horizontal dashed line 
in the same plot). This indicates that the electrical 
boundary conditions induce an enhancement in the po- 
larity of the film; this is unlike the vast majority of cases, 
where generally a suppression of P due to depolarizing 
effects is observed. The enhancement in P is due to the 
chemical bonding at the negatively polarized end of the 
film to the Pt surface. Such an effect was discussed for the 
centrosymmetric geometry in Ref. [t!; the present study of 
the fully relaxed capacitors in short circuit demonstrates 
that the effect persists in the polar structure. In fact 
an analysis of the local electrostatic potential shows that 
there is, instead of the usual depolarizing field (which 
would oppose the spontaneous P), a strong "polarizing" 
field; its magnitude is in excess of 200MV/m, and drives 
the film significantly more polar than the bulk. 

In such a "negative dead-layer" regime one would ex- 
pect thinner films to be even more polar. This is nicely 
confirmed by our results, shown in Fig. [51 where the po- 
larization enhancement attains values as large as 35% 
in the 2-cell thick capacitor; also the tetragonal ratio 
and the internal electric fields (not shown) steadily in- 
crease with decreasing thickness as expected from the 
above qualitative arguments. 
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FIG. 8: (Color online) Calculated values of the electric dis- 
placement field (or surface density of free charge stored on the 
plates) for short-circuited BaTiOa capacitors (black solid line 
and star symbols). The calculated bulk spontaneous polar- 
ization of BaTiOs at the SrTiOa in-plane lattice parameter is 
shown as a horizontal dashed line. Red squares are the values 
of the macroscopic polarization Pmac of the 5-unit cell and 
7-unit cell capacitors, as defined in Eq. (|27p . 



3. Layer polarizations and macroscopic polarization 

While it is not a necessary step to computing the 
macroscopic P, it is nonetheless interesting to push fur- 
ther the analogy to insulating superlattices, and compute 
the layer polarizations (LP) as defined in Ref. [l^. This 
involves grouping the Wannier centers of Fig. [7] into the 
separate clusters corresponding to the individual oxide 
layers, which are obvious from the plot, and computing 
the individual dipole moment pj per surface unit (see 
methods section). We plot in Fig. [5] the results for the 
Wannier-based layer polarizations. In all cases the dipole 
moment of the first BaO layer is about three times larger 
than the LP values in the rest of the films. This qual- 
itatively reflects the strong structural distortion, due to 
the chemical interaction with the Pt surface, which was 
discussed in the previous subsection. Otherwise, the LPs 
display qualitative features which are remarkably simi- 
lar to bulk BaTiOa. In particular, the LP of both layer 
types (BaO and Ti02) have the same positive sign (con- 
sistent with the positive displacement of all cations with 
respect to the O sublattice) , with the LPs of the BaO lay- 
ers systematically larger than the Ti02 values (approx- 
imately by a factor of 1.6). Note that in all cases the 
LPs are larger than the corresponding bulk values, and 
uniformly increase as the film becomes thinner. This con- 
firms that the enhanced structural distortions discussed 
in the previous section correspond indeed to an enhanced 
polarization of the films. 

It is interesting to note that the LPs converge rather 
quickly to a uniform bulk-like sawtooth pattern two or 
three oxide layers away from the interface, which indi- 
cates that the perturbations induced by the electrode are 
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FIG. 9: (Color online) Calculated Wannier-based layer po- 
larizations for the capacitors discussed in the text. Even- 
numbered layers are Ti02; odd layers are BaO. Bulk val- 
ues are reported for comparison as filled circles connected by 
dashed lines. 

rather local. This contrasts sharply with the picture pro- 
posed in a recent workii for the same system. We defer 
a detailed discussion of this issue to Sec. IVl 

The local LP values in the middle of the ferroelectric, 
PBaO and pTi02 J together with the average out-of-plane 
strain Cavg inferred from the data of Fig. [5l provide an 
accurate estimate of the macroscopic P inside the film as 

^'film = (pBaO + PTiOa) • (27) 

Cavg 

We plot Pfilm of the 5- and 7-unit cell capacitors as 
square symbols in Fig. [SJ here PfUm can be directly com- 
pared to the calculated values of the electric displace- 
ment D. First, the values of D and Pfiim are very close, 
as can be expected from their relationship (in SI units) 
D = eoffiim + -Pfilm: Indeed, eoE << P in typical ferro- 
electrics. Next, the fact that D is slightly larger than P 
is consistent with the electric field's being coUinear with 
P in the capacitors considered here, i.e., the interface in- 
duces a polarizing effect instead of a depolarizing one, as 
we already discussed extensively. 

Note that the electric field ffiim, as Pfiim, is the macro- 
scopically and planar-averaged value of E{x) inside the 
film and far from the interfaces. Therefore, we iden- 
tify £fiim and Pfilm as, respectively, the internal field and 
the polarization that are typically discussed in Landau- 
Ginzburg models of thin-film ferroelectrics. As will be- 
come clearer in the following sections, the relationship 
between ffiim, Pfiim and D is an intrinsic property of hulk 
BaTiOa, and does not depend on the interfaces, electrical 
boundary conditions or applied bias potential. This point 
is crucial to the development of our modeling strategies, 
and therefore we consider the above definitions of Pfiim 
and ffiim to be very convenient. Other author o^^'^^ de- 
fined P by averaging the dipoles over the whole volume 
of the film, including the interface region. Such a choice 
is less convenient for modeling, as i) it does not provide 
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a clear separation between bulk and interface effects; ii) 
it introduces a degree of arbitrariness, as the "bound" 
dipoles near the interface are strongly mixed with the 
metallic free carriers. 

Summarizing the above results, we have shown that, 
in thin-film capacitor configurations with Pt electrodes, 
BaO-terminated BaTiOs films display many features 
that are typical of a bulk-like crystal within the same 
mechanical boundary conditions (in-plane strain). How- 
ever, in addition to these similarities, there are also a few 
remarkable departures from the bulk behavior that are 
fully consistent with an interface-induced polarization en- 
hancement. The effect is stronger in thinner films, which 
is at odds with the common belief that realistic electrodes 
would systematically induce a polarization suppression 
due to imperfect screening. 

In the following we shall use the constrained-D method 
to understand the origin of this effect. In particular, 
we shall demonstrate that the "locality principle" is re- 
stored once the long-range electrostatic interactions are 
properly rationalized. In particular, we shall show that a 
simple model of this system with full ab-initio accuracy 
can be constructed in terms of the electrical properties 
of bulk BaTiOa and of the BaO-terminated Pt/BTO in- 
terface. 



C. Electrical equation of state 



ergy U. The results are plotted in Fig. [TOl Note the 
perfect match between the values of U calculated ab ini- 
tio (orange plus symbols) with the integrated potential 
(green curve) ; such a good match is a consequence of ac- 
curately compensating the Pulay stress with a fictitious 
constant negative pressure of tt = —2.61 GPa. The Pu- 
lay error is high (due to the relatively low plane- wave cut 
off of 40 Ry), and neglecting it would produce significant 
errors; with this simple correction, the numerical values 
are highly accurate. 

The relaxed ferroelectric ground state (e=0) is at 
rfo=0.364, which corresponds to P—SQAfiC/cm^, an en- 
ergy AC/=-28.6meV/cell, and c/a =1.061. Note that 
this is considerably larger than the value, c/a =1.038, 
in the strained centrosymmetric geometry. We can un- 
derstand the c/a ratio of the epitaxially constrained fer- 
roelectric phase as a result of two distinct contributing 
factors. One effect comes from the elastic properties 
(Poisson ratio) of the centrosymmetric crystal; straining 
BaTiOa to the SrTiOa lattice constant (—2.1%) produces 
a tetragonality of c/a =1.038 even in the absence of a po- 
lar distortion. The second effect, related to polarization- 
strain coupling, bring this value to c/a =1.061 once the 
unstable "soft" mode is condensed. For more details 
about the relationship between polarity and epitaxial 
strain we direct the reader to Ref. Hol. 



In order to model the electrical behavior of the 
Pt/BTO capacitors considered in this work, we now use 
the fixed- 13 method to sample the equation of state of the 
5-cell thick capacitor; this is the thinnest one in which 
Pcond is essentially zero in the middle of the oxide film, 
which ensures a high level of accuracy. For the scope of 
the present discussion, it is enough to restrict our investi- 
gation to a range of D that encompasses the equilibrium 
values calculated for different thicknesses in short circuit, 
i.e., 0.4e < d < 0.5e. Subsequently, we shall show how 
this information can be combined with the bulk equa- 
tion of state to predict the properties of capacitors of 
arbitrary thickness. We shall start by calculating the 
electrical equation of state of bulk BTO. 



1. Bulk BaTiOi 

As in the capacitor structures, we fix the in-plane lat- 
tice constant to ao = 7.276 a. u., the theoretical equilib- 
rium value for cubic SrTiOa, and we let the out-of-plane 
lattice parameter, as well as the internal coordinates of 
the five-atom tetragonal unit cell, relax as a function of 
the reduced electric displacement d. We use twelve val- 
ues of d, equally spaced between o?=0.0 and d — 0.55. 
We extract from the calculation the values of the po- 
tential and the c lattice parameter for every value of d; 
then we use splines to interpolate these values and we 
finally integrate the potential to recover the internal en- 



2. Capacitor structures 

We shall write the electrical equation of state of the 
5-unit cell capacitor as ez(d), the d-dependent reduced 
electric field. Since we already have one data point from 
the calculation in short circuit (d=0.478e, £=0), and we 
expect the potential to be a rather smooth function of 
d, it is likely that two additional points lying at the ex- 
tremes of the interval will be enough for constructing our 
model. Therefore, we repeat the calculation of the 5-cell 
capacitor twice with d set to 0.4e and 0.5e, respectively. 
(In both cases the ionic positions and out-of-plane lattice 
parameter are relaxed to the same convergence thresholds 
used in the zero-field cases.) The smoothness of the po- 
tential is confirmed by our results, plotted as black circles 
in Fig. [TTl which lie almost on a straight line. 

To account for the small curvature, we interpolate the 
points with a second-order polynomial expanded around 
the spontaneous reduced displacement do of the bulk. 
To establish notation, we do this first for a simple bulk 
crystal. Recalling that internal energies are related to 
the reduced electric fields by Eq. ([7|), e(d) = dU/dd, and 
anticipating an expansion of the internal energy up to 
third order in d — do, 

C/b(d) =A^o+ A\id - do) + ^(d - do)2 + ^(d- do)3, 

(28) 
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Reduced electric displacement d 



FIG. 10; (Color online) Internal energy, reduced electric field, and c/a ratio vs. reduced displacement field (in units of e) for 
coherently strained bulk BaTiOa. Ab-initio data are shown as symbols; dashed curves in the middle and right panels are spline 
interpolations; continuous curve in the left panel is the numerical integral of the spline in the middle panel. 



we expand e\i{d) as 

e^{d)=A\ + A\{d-dQ) + ^{d-dQf. (29) 

(Note that the expansion is carried out about the spon- 
taneous displacement do of the bulk, so that A}l vanishes 
by construction.) We then carry out similar expansions 
for each A^-cell capacitor, e.g., for Af=5, 

U,{d) = A'iHAf\d-doH^-{d^d,f + t^{d^dof 

(30) 
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FIG. 11: (Color online) Calculated values of reduced elec- 
tric field e as a function of d in capacitors of various thick- 
nesses (empty symbols). Solid lines correspond to the first- 
principles-derived model described in the text. Dashed line 
is the calculated bulk equation of state (from middle panel of 
Fig. nop , which corresponds to the potential drop across a sin- 
gle bulk unit cell. Star corresponds to the relaxed ferroelectric 
state of the strained bulk crystal. 



and 

4(5) 

e5(d) - Af + Af{d -d^) + ^{d- d^f. (31) 

The fitted bulk expansion parameters , and the inter- 
face parameters defined via 

A\ = 4f ) - hAl, (32) 

are reported in Table HI 

Focusing first on Eq. (|3ip for the 5-ccll-thick capacitor, 
the fitted £5(1^) is shovifn as the solid black line passing 
through the circles in Fig. [TT] Using this together with 
the bulk information encoded in Eq. ([^5)1 . we can then 
'predict the equations of state for thinner or thicker ca- 
pacitors according to the formula 

ejv(d) = (A^-5)eb(rf)+£5(d). (33) 

Setting TV to 2, 3 and 7 we obtain the colored solid curves 
in Fig. [TT] The intersection of each curve with the e = 
axis yields a well-defined value of d, which is the pre- 
dicted polarization state of a capacitor with thickness A^, 
and can be directly compared to the first-principles data 
already in hand. To that end, we take the d values from 
Fig. [5] and plot them as colored symbols on the e = axis 
of Fig. [TT] The agreement between the first-principles 
points and the model predictions is extraordinarily good 
(discrepancies are smaller than 0.1%). Surprisingly, this 
holds true even for the thinnest (2-cell) capacitor, where 
one would expect the estimation of d to be less accurate 
(see previous sections). Also, apart from purely techni- 
cal issues in defining c?, one might expect the properties 
of such a thin layer of oxide to depart somewhat from 
what is calculated in thicker capacitors. The accuracy of 
our model in this thickness regime is indeed encouraging, 
and indicates that our methods for accessing the inter- 
facial electrical properties are able to predict, with full 
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FIG. 12: (Color online) Calculated energies, relative to corre- 
sponding paraelectric state, for BaTiOa (BTO) capacitors of 
various thicknesses and polarization states (symbols). Solid 
curves correspond to the model discussed in the text. Verti- 
cal dashed line indicates the spontaneous polarization of bulk 
BaTiOa. 

first-principles accuracy, the behavior of a wide range of 
systems. 

To further confirm the internal consistency of our 
model, we perform a similar analysis for the energetics 
of the capacitors. Let Uj\r{d) be the difference in inter- 
nal energy, for a given thickness N, between the state 
at specified d and the paraelectric structure at d=0. We 
plot in Fig. [13 the three values of (/^{d) (black circles) 
that we extracted for the 5-cell capacitor from the same 
calculations described above. Expanding in d ~ do ac- 
cording to Eq. we note that all the coefficients have 
already been determined from Eq. (|3ip except for the 

(5) 

arbitrary constant of integration Aq . Adjusting this 
one free parameter, we find an excellent match of the fit 
(black curve) with all of the data (black circles). As was 
done for the potential, we then predict the energy UM{d) 
for other N just by adding or subtracting bulk units, 

UN{d) = {N- 5)Ui,{d) + U5{d)- (34) 

Again, all of the needed bulk coefficients were already 
determined from Eq. except for the constant term 
= AU. The results for N^2, 3 and 7 are plotted as 
the colored curves in Fig.[T21 Again, the minima of all the 
C/7v(rf) curves match very well the points explicitly cal- 
culated in the short-circuit first principles calculations at 
various thicknesses. The maximum discrepancy is about 
1 meV, which is comparable to the numerical noise in the 
total energy values introduced by the discreteness of the 
plane- wave basis set during variable-cell structural relax- 
ations. 

Two important details are apparent from Fig. [11] and 
Fig. [m First, all curves in Fig. [TT] have a common inter- 
section at d = do; this is related to the fact that at d — do 
the internal electric field in the bulk vanishes: A\ = (we 
shall come back to this point in Section fill E[) . Second, 
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TABLE I: Values, in atomic units, of the expansion coeffi- 
cients used to model the bulk and interface contributions to 
the electrical equation of state, Eqs. (|28l29p and p6l37|l , and 
to the elastic equation of state, Eqs. (140 1) and (|43|) . 



the relaxed internal energies of the capacitors considered 
here are about one order of magnitude larger than the 
depth of the bulk double well (see Fig. [TH). This indi- 
cates that the chemical bonding mechanism discussed in 
Ref. |7| has a substantial impact on the energetics, which 
explains the strong tendency of the capacitors towards a 
superpolar state. 

D. Interfacial dielectric and piezoelectric response 

Our goal now is to show how several useful interface- 
specific observables can be directly linked to the electric 
equations of state (EOS) discussed above. In addition to 
the purely electrical variables, we shall further extend our 
model by addressing also the elastic (i.e., piezoelectric) 
EOS of both bulk and electrode interface. 



1. Dielectric response 

The polynomial expansion of the dielectric response 
(electrical EOS) was essentially already determined in 
Sec. lIII C"2l Using the = 5 as our model structure and 
using Eq. (|33p . we single out the interface contribution 
by defining the interface EOS to be that of a hypothetical 
"zero-thickness capacitor" , 

ei(d) = eo{d) = £z{d) - 5ebuik(d), (35) 

with a similar relation for the internal energy. The inter- 
face potential and energy are then expanded as in analogy 
to Eqs. and Eqs. ([3DHSI1) as 

A A^ 

UM^Y.^i'^-^oT' (36) 

n=0 

s-iid)^Y.^id-dor. (37) 

Tl = 

The coefficients Ajj are determined once and for all from 
a pair of calculations on the bulk and on the 5-cell capac- 
itor superlattice using Eq. ([5^ . The resulting bulk and 
interface coefficients are reported in Table ID 
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The physical interpretation of the zero-order coefRcient 
Aq is immediate, as it represents the interfacial contribu- 
tion to the energy of the capacitor when a reduced dis- 
placement do is induced in the film (the energy zero is 
set to that of the paraelectric d = state). The first- 
order coefficient is also physically transparent, as it cor- 
responds to the interfacial potential drop at d = do- As 
we mentioned above, the bulk ground state at has zero 
internal field, so an applied external bias of A\ will al- 
ways induce the same (bulk-like) polarization, regardless 
of the thickness N of the film. 

In order to interpret the higher-order coefficients, we 
first derive an expression for the inverse capacitance that 
is valid for bulk, interface and the full capacitor structure. 



C^'id) 



de^jd) 
dd 



Ai 



Af{d- 



do) 



(38) 



where X=b,I,(A^). Therefore, A^ is the inverse capaci- 
tance at d — do- In the bulk case we can directly link this 
coefficient'-- to the static (free-stress) dielectric constant 
of ferroelectric BaTiOa through Eq. pT|) . 



^33 



47rc 
SA^- 



(39) 



Using the values reported in Tab. [J we obtain €33 — 37. 

In the capacitor case the physical meaning of C^~J^^ (d) 
is obvious; note that this is an inverse capacitance per 
surface unit cell and must be multiplied by S to obtain 
an inverse capacitance density. Finally, Cj~^(d) is the 
inverse interface capacitance that was discussed, for in- 
stance, in Ref. 0- (Note that here Cj is the combined 
effect of both interfaces, unlike the C/ defined in Ref. 7. 
We shall present a general strategy for separating this 
quantity into two individual contributions later, in Sec- 
tion |TVl) By using the values of Table U we obtain an 
interfacial capacitance density Ci{do)/S = 0.44 F/m^; 
this value, due to the dielectric non-linearity contained 
in the third-order coefficient A\, is reduced by half near 
the right end of the d interval considered here {d ~ 0.55). 
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FIG. 13: Calculated interfacial distance c^(d) of Eq. (|42)) for 
capacitors of various thicknesses and polarization states (sym- 
bols); solid curve is a quadratic fit. 



With the calculated values, we obtain ^33 = 30 pm/V. 
(Note that in our calculations we kept the in-plane lat- 
tice parameters fixed; this value might change upon full 
relaxation of the crystal.) 

Next, we define the interface contribution c^{d) as 



c^{d) = L]y{d) - iVcbuik(rf) - nptS 



pt- 



(42) 



where Lis[{d) is the total thickness of the relaxed capaci- 
tor for a given value of d, and N and npt are the number 
of oxide and metal bulk cells in the capacitor superlattice. 
Note that Spt, the relaxed interlayer distance in bulk Pt 
(strained in plane to the same lattice parameter), is inde- 
pendent of d. Using this formula, we extracted c^{d) for 
all the capacitor structures considered so far, and plotted 
the values in Fig. [T31 In the same figure the solid line is 
a quadratic fit using the formula 



c'{d) 



ciid- 



do) + f (d- 



do) 



(43) 



2. Piezoelectric response 

In order to describe piezoelectric effects, we consider 
the bulk lattice parameter Cb and the interfacial distance 

as functions of d. In the same spirit as before, we first 
perform a quadratic fit of Cb as 



Cb(d) -c^ + c^(d-do) + ^(d-rfo)' 



(40) 



Here Cq is the equilibrium lattice parameter of the epi- 
taxially strained tetragonal state, and c\ is related to the 
piezoelectric constant by Eq. ((TO)) . 



dcb 

"33 = -7^ 

de 



d—dn 



dCh 
dd 



de\-i 
dd) 



d—dn 



The zeroth-order coefficient Cg has the meaning of a com- 
bined effective interface distance for both top and bottom 
bottom electrodes. (It may look larger than expected be- 
cause there are actually iV-t- 1/2 oxide cells in our capaci- 
tors, not N, and because the Pt/BTO interface distances 
are typically somewhat larger than the bulk interlayer 
spacings of either Pt or BTO.) The linear coefficient, in 
analogy to the bulk case, is related to the electrode con- 
tribution to the piezoresponse of the capacitor. All the 
coefficients defined in the text are reported in Tab. HI 

In summary, our simple model is able to predict, within 
the numerical accuracy of a full first-principles calcula- 
tion, the energy, polarization, dielectric and piezoelectric 
response of a Pt/BTO/Pt capacitor of arbitrary thick- 
ness (from two unit cells to infinity). 
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E. Band lineup and ferroelectricity 

The simple model derived above allows us to inter- 
pret the physics behind the polar enhancement from yet 
another point of view. If we go back to Fig. [Til it is ap- 
parent that the solid lines converge to the same point 
at a value of d corresponding to the spontaneous po- 
larization of bulk BaTiOa {do = 0.364). The reason is 
that the potential is zero in BaTiOs at that value of d, 
because the bulk crystal is at electrostatic equilibrium; 
hence, adding or subtracting bulk units does not change 
the value of £Ar(do)- Therefore, e/(rfo) — £^(0^0), inde- 
pendent of N, is an intrinsic interface property related 
to the band lineup between the metal electrode and the 
fully polarized ferroelectric film, 

£/(do) -4 =0+ -0-, (44) 

where cj)^ are the respective Schottky barrier heights 
at the positively and negatively polarized ends of the 
film. Since the symmetry is broken upon ferroelectric 
off-centering, these two values generally differ. 

Based on typical assumptions of phenomenological the- 
ories and on a large body of experimental data, one would 
expect a positive interface potential ei(da) > 0; in short 
circuit this would yield a depolarizing field that opposes 
the polarity of the film (recall that e = —V). This means 
that generally one needs to apply a positive bias in or- 
der to reach (and sustain) bulk values of d in a thin 
film; when the bias is switched off, the polarization is 
either reduced or relaxes to zero by transitioning to a 
multidomain state. "^^ The Pt/BTO/Pt system analyzed 
in this work displays a rather different behavior, in that 
£i{do) ~ —0.8 V is negative. This is a signature of the 
strong polarizing field discussed earlier in the context of 
our calculations in short circuit. 

Note that in our calculations so far we never extracted 
(j)'^ and separately, since in the case of compositionally 
symmetric capacitors only their difference matters for the 
electrical properties of the device. Analogously, we con- 
sidered only a total interfacial distance c/ that takes into 
account both the top and bottom ends of the film. In the 
following sections wc shall use the techniques described 
in Section HID 21 to single out the potential lineup and 
distances of either interface individually as a function of 
d. We shall demonstrate, as an example, how this in- 
formation allows one to accurately model the electrical 
behavior of truly asymmetric devices starting from calcu- 
lations performed on systems having a centrosymmctric 
paraelectric geometry. 

IV. RESULTS: Au/BaZrOs/Au CAPACITORS 

A. Motivation 

As a model system for studying the electrical behavior 
of asymmetric capacitors we choose ferroelectric devices 



with Au as the electrode and BaZrOa (BZO) as the ac- 
tive film. This choice is motivated by issues of compu- 
tational practicality. As mentioned in the introductory 
sections, an essential prerequisite for defining and con- 
trolling the polarization in a metal/insulator heterostruc- 
ture is its insulating character along the direction perpen- 
dicular to the interface. This can only be true if there 
are non-vanishing Schottky barriers at both interfaces, 
and if those barriers are preserved upon ferroelectric off- 
centering. The wider band gap of BZO makes this prop- 
erty much easier to satisfy than with more conventional 
ferroelectric oxides such as PbTiOs and BaTiOs. Then, 
given the larger lattice parameter of BZO, we decided to 
use Au as the electrode instead of Pt in order to avoid 
unrealistically large strains in the electrode slab. (Pt is 
more popular in the ferroelectrics community because it 
matches better the lattice parameter of many Ti-based 
perovskites.) This combination of materials yields large 
Schottky barriers for both BO2- and AO- terminated in- 
terfaces and is therefore ideally suited to the scope of the 
present study. An interesting aspect of this study is that 
we are able to compare the electrical behavior of BO2- 
and AO-terminated interfaces. 

One important issue is that BZO is not ferroelectric. 
Experimentally it has a stable cubic structure down to 
very low temperatures, while theoretically there have 
been several reports of zone-boundary instabilities asso- 
ciated with rotations of the oxygen octahedra.— In or- 
der to induce a polar instability in BZO, we set the in- 
plane lattice parameter to a fixed value of 7.60 a. u. (a 
compressive strain of —3.0% with respect to the theo- 
retical equilibrium lattice parameter of 7.38 a. u. of the 
cubic structure).'^? Note that we did not check for possi- 
ble competing non-polar states, as such an analysis would 
require doubling the size of the simulation cell, substan- 
tially raising the computational cost. For this reason, 
our setup should not be understood as a direct prediction 
of ferroelectric behavior in epitaxially strained BaZrOs. 
Rather, we intend it primarily as a tractable computa- 
tional model, which we expect may be representative of 
the behavior of a typical perovskite with a tetragonal fer- 
roelectric ground state (e.g., PbTiOa or BaTiOs at room 
temperature). 



B. Computational model 

Our model heterostructures consist of (OOl)-oriented 
BZO films with symmetrical Zr02 or BaO terminations 
and a thickness of 8.5 unit cells, interfaced with a metal 
electrode slab of 11 Au monolayers. Our first goal is 
to study the full equations of state of these structurally 
symmetric capacitors, to establish the similarities and 
the differences arising from the dissimilar (Zr02/Au vs. 
BaO/Au) bonding configurations. Next we extract the 
interface-specific information and use it to predict the full 
equation of state of an asymmetric configuration (with 
BaO and Zr02 terminations at opposite ends). Then, we 
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FIG. 14: (Color online) Internal energy, reduced electric field, and c/a ratio vs. reduced displacement field (in units of e) for 
coherently strained bulk BaZrOa. Ab-initio data are shown as symbols; dashed curves in the middle and right panels are spline 
interpolations; continuous curve in the left panel is the numerical integral of the spline in the middle panel. 



verify that our procedure yields the desired result by com- 
paring this prediction with the explicitly computed U{d) 
curve for an asymmetric capacitor having 8 unit cells of 
BZO and 12 layers of Au. Before going into details about 
the capacitor structures, however, we first briefly summa- 
rize the electrical properties of bulk BaZrOa within the 
symmetry and mechanical constraints described above. 

C. Bulk BaZrOa 

As in the capacitor structures, we fix the in-plane lat- 
tice constant to oq = 7.60 a.u. and let the out-of-plane 
lattice parameter, as well as the internal coordinates of 
the five-atom tetragonal unit cell, relax as a function of 
the reduced electric displacement d. We use seven evenly 
spaced values of d ranging between d—O.Q and (i=0.276e. 
We extract from the calculation the values of the reduced 
field £ and the lattice parameter c for each value of d, use 
splines to interpolate these values, and finally integrate 
e{d) with respect to d to recover the internal energy U{d). 
The results are reported in Fig. [141 As for the BaTiOs 
case, the match between the values of U calculated di- 
rectly (red 'plus' symbols) with those obtained by inte- 
grating £{d) (black curve) is excellent. Both the sponta- 
neous polarization and the double- well potential depth 
AU are significantly smaller than in BaTiOa. The re- 
laxed ferroelectric ground state is at c?=0.2076, which cor- 
responds to P=20.5/iC/cm2, AC/=-3.05meV/cell, and 
c/a =1.052 (compared to c/a =1.0475 in the strained 
centrosymmetric geometry). 

D. Schottky barriers 

In Fig. [15] we plot, as a function of the reduced displace- 
ment d, the Schottky barrier heights (SBH) extracted 



from the symmetric BaO- and Zr02-terminated capaci- 
tors using the techniques of Sec. Ill D 21 All values lie be- 
tween about — 1.2eV and — 1.7eV. Considering that the 
calculated LDA gap for centrosymmetric bulk BaZrOs is 
3.12 eV, this indicates that in all cases the Fermi level 
of Au lies close to mid-gap, and that our BZO/Au ca- 
pacitors are thus free from Schottky-breakdown issues 
Both curves shown in Fig. [T5| have considerable curva- 
ture, although with dissimilar features. In the case of 
Zr02/Au this curvature is found to be dominated by 
the bulk AVF(d): If we remove such a dependence from 
0ZrO {d) , we obtain a roughly linear function (red dashed 
curve, which is related modulo a constant shift to the in- 
terfacial step in the electrostatic potential). By contrast, 
removal of the bulk contribution /Sy-p{d) (black dashed 
curve) does not restore linearity in the BaO/Au case. 



Incidentally, we note that our method may be of gen- 
eral use for the computation of Schottky barriers^^*^ (in- 
dependently from their relationship to dielectric proper- 
ties) , which are extremely important for many technolog- 
ical applications. Recall that the SBH has a truly unique 
definition only when the macroscopic electric field van- 
ishes in the insulator, since then one can identify the 
valence and conduction band edges precisely. For non- 
centrosymmetric insulators such as wurtzite oxides or 
spontaneously polarized ferroelectrics, such a condition is 
not easily obtained in ordinary first-principles supercell 
calculations. However, our approach makes it extremely 
easy to do such calculations. For example, we have indi- 
cated with large 'star' symbols in Fig.[T5|the "physical" 
values of the SBH, i.e., those corresponding to the spon- 
taneously polarized ferroelectric film in the absence of 
any internal field. 
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FIG. 15: (Color online) Solid curves: Interfacial p-type Schot- 
tky barrier values extracted from the symmetric capacitor 
configurations for BaO-terminated (black circles) and ZrO- 
terminated (red squares) films. Dashed curves: Same, except 
that the bulk dependence of the VBM on the electric dis- 
placement has been subtracted out. The sign of d follows the 
convention that the electrode lies at « > 0. 



E. Interfacial equation of state: Zr02/Au 

The discussion in the previous section suggests that 
there might be important qualitative differences between 
the behavior of BaO/Au and Zr02/Au interfaces. How- 
ever, it is not immediately obvious how to interpret the 
4>{d) curves directly, as their relationship to the physical 
electrical response of the capacitor contains some aspects 
of arbitrariness. Such arbitrariness does, of course, cancel 
out when the final equation of state of the entire capaci- 
tor is constructed. Therefore, in order to obtain quanti- 
ties that have a direct physical meaning, we proceed by 
combining the above 4>{d) curves in pairs as appropriate 
for the capacitor structures of interest. There are four 
such structures that we denote as 'AB,' where A and B 
are variables that specify, for the bottom and top inter- 
faces respectively, whether the interface is BaO/Au or 
Zr02/Au. Then the interface contribution to the equa- 
tion of state of the specified capacitor structure is 

£i.AB(cf) = -0A(-d) + 0B(d), (45) 

in terms of which the the equation of state of the entire 
TV-cell capacitor is then eN,AB{d) — ei,AB{d) + Neh{d). 

The four resulting functions ei.AB(c?) are plotted in 
Fig. [ini The most striking feature is the almost per- 
fect linearity of the symmetric Zr02/Zr02 configuration. 
This means that the description of the interfacial equa- 
tion of state in terms of a constant interfacial capaci- 
tance (i.e., replacing the interfaces by a layer of linear 
dielectric in series with a bulk-like BaZrOs film) is ap- 
propriate in this case. The slope of the ei{d) curve yields 
a combined capacitance density of Cj/ S=1.73 F/m^ for 
both interfaces, so that each interface is associated with 
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FIG. 16: (Color online) Interfacial equations of state, re- 
constructed from the Schottky barrier values of Fig. 1151 for 
a symmetric BaO- (black circles) or Zr02-terminated (red 
squares) capacitor. The blue upward-oriented triangles re- 
fer to an asymmetric arrangement with the BaO termination 
on top; the green downward-oriented triangles have the BaO 
termination at the bottom. 

a capacitance density of 3.46 F/m^. This value is re- 
markably high when compared to the typical range of 
~0.4-0.6 F/m^ calculated^i2§' for oxide electrodes such as 
SrRuOs. This result corroborates the ideas proposed in 
Ref. Ill that weak electrode-oxide bonding is beneficial to 
the ferroelectric properties of a capacitor. Here we in- 
deed find that a chemically inert electrode material such 
as Au yields excellent screening and only a marginal per- 
turbation to the polar response of the film. Using the 
formalism developed in Ref. 0, we find a "critical thick- 
ness for ferroelectricity" iVcrit = 3 for both symmetric 
geometries (Zr02/Zr02 and BaO/BaO). 



F. BaO/Au and bonding properties 

Interestingly, the BaO/BaO curve is almost exactly 
overlapping with the Zr02/Zr02 one in the interval 
— 0.15e < d < 0.15e, while a strong departure from the 
linear regime occurs for values of d lying outside this in- 
terval. A significant non-linearity was indeed expected 
from the (f>{d) curves in Fig. [151 the fact that this non- 
linearity cancels for — 0.15e < d < 0.15e and yields a 
quasi-linear behavior is probably coincidental. The non- 
linearity of the BaO/Au interface emerges most clearly 
in the case of an asymmetric capacitor (green and blue 
curves in Fig. 1161 which are correctly related by a mirror 
symmetry operation). 

To trace the origin of the qualitative difference between 
Zr02/Au and BaO/Au interfaces (linear vs. non- linear 
behavior), it is useful to follow the evolution of the Au-0 
bond length as a function of electric displacement, plot- 
ted in Fig. [T71 While the bond length varies only weakly 
and follows a linear trend for Zr02/Au, it covers a much 
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FIG. 17: (Color online) Interfacial Au-O bond distance as a 
function of electric displacement field. The sign of d follows 
the convention that the electrode lies at z > 0. 

wider range of distances (2.2 to 3.1 A) and displays a 
strong nonlinearity for BaO/Au. We interpret the latter 
behavior as indicative of breaking and reforming of the 
Au-O bond upon polarity switching. Clearly, the break- 
ing of a bond is a highly non- linear event, helping to 
explain the calculated features of the electrical response. 
This picture agrees fully with the arguments of Ref. 7, 
where the bond stability (instability) was correlated with 
the suppression (enhancement) of the tendency to polar 
instability of the capacitor structure. The present results 
corroborate these ideas and provide further evidence for 
the strong correlation between interfacial chemistry and 
electrical response. 

Interestingly, while in the case of Pt/BaTiOs/Pt the 
interface bonding mechanism strongly enhances ferro- 
electricity, in the present case of Au/BaZrOa/Au the 
overall effect is a slight suppression. We shall briefly dis- 
cuss the origin of this dissimilar behavior in the following 
subsection. 



G. Enhancement or suppression? 

Why do certain AO-terminated interfaces (e.g., 
BaTiOa/Pt, BaTiOa/Au) enhance the ferroelectric in- 
stability of the film, while others (especially PbTiOa/Pt 
but also, to a smaller extent, the BaZrOa/Au one dis- 
cussed here) suppress it instead? While a definitive an- 
swer is not yet available, some qualitative trends can be 
explained in terms of the frustrated bonding-environment 
model of Ref. 'tI. According to this model, a flat AO layer 
in contact with the electrode produces a competition be- 
tween the A-metal repulsion and the 0-metal attraction. 
The buckling of the AO layer caused by a bulk ferroelec- 
tric distortion of one sign or the other shifts the balance, 
causing the bonding or the repulsive force to prevail. In 
fact, even in the centrosymmetric capacitor geometry, the 



interface AO layer is not flat, but exhibits a certain de- 
gree of buckling (with the A cation typically displacing 
toward the oxide film) due to the broken-symmetry en- 
vironment. One can therefore expect some difference in 
behavior between perovskite AO-terminated films that 
show different degrees of "natural buckling" at the sur- 
face of their cubic reference phase. We computed the 
values of the AO rumpling of the free PbTiOa, BaZrOa 
and BaTiOa surfaces, finding values of 0.136, 0.121, and 
0.022 A, respectively. Indeed, these results indicate that 
the film with the much flatter surface (BaTiOa) displays 
a strong enhancement of the polar instability, while those 
that are significantly buckled (PbTiOa and BaZrOa) do 
not. While this is only a rough indication, and other fac- 
tors are most likely at work, it suggests a correlation that 
may help to explain our detailed numerical results. 

H. From symmetric to asymmetric 

We claimed earlier that it should be possible to use 
the interface equations of state extracted from calcula- 
tions on symmetric capacitors to predict the equation 
of state for the asymmetric case. We demonstrate this 
now. Our "asymmetric" geometry is comprised of an 8- 
unit-cell BaZrOa film and a 12-layer Au slab, where the 
bottom and top interfaces (relative to the oxide) are of 
Au-BaO and Zr02-Au type respectively. Thus, a positive 
d corresponds to the polarization pointing towards the 
Zr02-Au interface. The interface-specific contribution to 
the reduced electric field, £/(c?) — —(j)L{d)+4>R{d), is plot- 
ted as the green curve in Fig. [161 Note that the curve is 
linear for d > 0, while it shows a significant non-linearity 
for d < 0, where the Pt-0 bond breaks, in agreement 
with the discussion of the previous sections. We now use 
this function to reconstruct the e(d) curve of the whole 
capacitor by adding an appropriate number of bulk units 
as in Eq. (ffg)) . 

£(d) =e/(d) + 7Vebuik(d), (46) 

with N — 7.5."^^ Then we numerically integrate e{d) to 
obtain the U (d) energy curve, which is plotted as the 
dashed green curve in Fig. [TSl 

In the same figure we plot two other curves correspond- 
ing to the symmetric geometries, which were obtained 
from the symmetric curves in Fig. I16p in the correspond- 
ing way (with iV = 8). For all three curves we also plot, 
as symbols, the values of U extracted directly from the 
first-principles calculations. The match is almost per- 
fect for the symmetric cases, as expected since the po- 
tentials and energies were taken from the same calcu- 
lation; their agreement is just a test of internal consis- 
tency. (Incidentally, note the close agreement between 
red and black curves in Fig. [T51 especially in the cen- 
tral region, which is inherited from the similarity of the 
corresponding red and black curves in Fig. [111) The 
acid test, however, concerns the asymmetric structure: 
the points were extracted from a direct calculation on 
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FIG. 18: (Color online) Electric equation of state U{D) for 
three capacitor configurations discussed in the text. Curves 
correspond to integrating Eq. p9p [or Eq. (|46p ]. where 0L,R(rf) 
were extracted from the symmetric capacitor calculations 
(solid curves in Fig. I15|l . and ebuik(d) was extracted from 
the bulk calculation (middle panel of Fig. I14p . Points corre- 
spond to the U{d) values explicitly calculated for symmetric 
BaO:BaO (black squares), symmetric Zr02:Zr02 (red circles) 
and asymmetric BaO:Zr02 (green diamonds) capacitors. 

the asymmetric structure, while the curve was inferred 
from the data on symmetric capacitors using our model 
of Eq. p^ . The match is excellent, with less than a 1% 
discrepancy in the spontaneous polarization at the mini- 
mum in Fig. llSi corresponding to short-circuit boundary 
conditions. (The predicted and calculated D values are 
—22.1 and — 22.2/iC/cm^ respectively.) Note that the 
spontaneous D is significantly enhanced (~ 8%) com- 
pared to the bulk value. This is a consequence of the 
nonlinearity discussed above, which induces a positive £/ 
for d < (see green curve in Fig. [TB)) . The prediction 
and the actual calculation also nicely agree regarding the 
absence of a secondary minimum at positive d. Inter- 
estingly, both strategies would yield such a minimum if 
the capacitor were just one unit cell thicker; this further 
confirms the accuracy of our model. 

While all these physical aspects compare very favor- 
ably, Fig. [T51 shows also some discrepancies in the actual 
values of the internal energy U, in particular concerning 
the depth of the energy minimum (predicted and cal- 
culated AU values of -38.4 and -41.4meV respectively). 
We shall briefiy discuss this discrepancy in the following 
Section. 



I. Accuracy issues 

We demonstrated in Section HTTl that the "locality prin- 
ciple" established in Section HIl holds very accurately, al- 
lowing one to predict the electrical properties and ener- 
getics of capacitors of varying thickness with excellent 
fidelity. In this respect the discrepancy in the energy 



curves of Fig. [18] is somewhat surprising, and it is worth 
discussing its origin to make sure that all aspects of the 
method are under control. 

We find that the cause of the discrepancy is rooted 
in the treatment of the Au electrode in the simulations, 
rather than in numerical (or formal) errors. In the sym- 
metric calculation, we used an 11-layer Au slab, and from 
this data we constructed the green curve in Fig. [THl On 
the other hand, we had to choose an even (12-layer) Au 
slab in the asymmetric calculation, which yielded the 
green diamond symbols. It is well known that quantum 
size effects associated with Fermi-surface nesting can per- 
sist to substantial thicknesses in thin metal films4§- Thus, 
it is not unreasonable to expect the surface of a 12-layer 
slab to behave slightly differently from that of an 11- 
layer slab. To prove this point, we calculated the work 
function of two free-standing slabs of 11 and 12 layers, 
and we found a difference of about 60 mV - enough to 
produce non-negligible discrepancies in the capacitor cal- 
culations. The use of a finite electronic temperature (to 
accelerate convergence with respect to fc-point sampling) 
helps, in that it makes the one-particle density matrix of 
the metal short-ranged in space. However, further explo- 
ration of these issues falls outside the main scope of the 
present work, and we satisfy ourselves with advising the 
reader of this issue so that it can be kept in mind when 
performing future calculations. 

V. DISCUSSION 

In the following, we discuss the implications of our 
work by comparing our techniques and results to the rel- 
evant literature. 



A. Locality 

One of the crucial aspects of our work concerns the use 
of the "locality principle" in the simulation of capacitor 
structures. As we mentioned in the introduction, at least 
two recent theoretical works have reported phenomena 
which do not appear consistent with this assumption. 
We shall briefiy discuss them here. 

First, the authors of Ref. [ill reported, for BaO- 
terminated Pt/BaTiOs/Pt capacitors (structurally anal- 
ogous to those considered in Section HII)) . a "ferrielectric" 
layer-polarization (LP) pattern, with profound quali- 
tative deviations from the bulk pattern, affecting the 
whole volume of the oxide film. This in principle im- 
plies a sharp, qualitative, deviation from the "locality 
principle" mentioned above. While we cannot provide 
a definitive explanation for the origin of the disagree- 
ment with our results, we believe it might lie in the 
subtleties involved in the LP construction for an over- 
all metallic system such as the ferroelectric capacitors 
under consideration. In contrast with our work, where 
localization of the Wannier states is imposed in one di- 
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mension separately for each fc-point (as in the origi- 
nal LP formulation in Ref. I19f ). the authors of Ref. [Ill 
used fully three-dimensional maximally-localized Wan- 
nier functions. Furthermore, they treated the metallic 
states rather differently than here, in that they adopted 
a preliminary disentanglemenli^ procedure before local- 
izing the states by means of the Marzari-Vanderbilt algo- 
rithm.^- Finally, the authors of Ref. [ill might have taken 
different prescriptions for the assignment of the Wannier 
functions to the individual oxide layers, possibly intro- 
ducing the LP equivalent of the quantum of polarization 
in their reported values. Regardless of which of the above 
factors may be responsible for the discrepancy, the effect 
proposed in Ref. [TH appears likely to be a consequence 
of the details of the Wannier localization/grouping pro- 
cedure, and we therefore suggest that its physical signif- 
icance should be judged with some caution. 

Second, the authors of Ref. [l^ considered 
SrRuOa/KNbOs/SrRuOa capacitors and reported 
an interface-induced disruption of the ferroelectric soft 
mode of the film, with the appearance of a head-to-head 
"interface domain wall" located three unit cells away 
from the electrode. This was interpreted as an effect 
of the strong bonding at the interface, which would 
"clamp" the interface dipoles to a fixed value; this 
constraint would then couple with the ferroelectric 
instability of the film, producing the calculated inhomo- 
geneous polarization pattern. The spatial variation in 
P is strongly asymmetric, and takes place over approxi- 
mately 6 unit cells at the positively polarized end of the 
KNbOs filmiiS This contrasts with our results, where the 
interface-induced distortions are extremely local and heal 
completely within the first perovskite unit cell adjacent 
to the interface. Even if the metal-insulator interactions 
are somehow stronger for the SrRuOa/KNbOs system 
than for our cases, this should just be reflected in a 
stronger functional dependence in the interface equation 
of state, modifying the strength of the depolarizing field. 
One should still expect a uniform polarization deep in 
the insulator, unlike the inhomogeneous polar ground 
state found by these authors. Therefore, it is difhcult to 
understand their findings unless one of the fundamental 
prerequisites for the formalism developed in this work 
might have broken down. For example, Junquera and 
Ghoseai have emphasized the dangers of pathological 
band alignments which, as an artifact of the band-gap 
problem of density-functional theory, may lead to charge 
spillage into the perovskite at certain perovskite-metal 
interfaces. We suggest that such a possibility should 
be investigated for the SrRuOs/KNbOa interfaces 
considered in Ref. [Tol . 



B. Relationship to Landau theory 

Our approach has many points of contact with ear- 
lier Landau-theory models of depolarization in thin-film 
ferroelectrics. However, there are also some notable dif- 



ferences that we shall emphasize in the following, in order 
to avoid confusion or misunderstanding. 

First, we note that our strategy is different in spirit 
from what was done, for example, in Ref. There, 
the authors fitted the parameters of a Landau-like ex- 
pansion to the calculated first-principles values of the 
depolarizing field in short-circuited capacitors. In con- 
trast, in our approach the first-principles engine works 
as a stand-alone tool that yields the ground state en- 
ergy and structure as a function of a well-defined elec- 
trical variable (within a given set of specified mechan- 
ical/symmetry constraints and thermodynamic ensem- 
ble). These data can then be fitted a posteriori to a 
polynomial, thus obtaining an expression that bears a 
close resemblance to Landau-theory expansions, but the 
latter is not a necessary step. 

Another difference concerns our use of an interfacial 
capacitance (or, equivalently, of an effective screening 
length) that embodies all the physical ingredients con- 
tributing to the electrostatics. Gerra et air^ on the 
contrary, made a distinction between purely electrostatic 
screening and short-range chemical bonding effects, and 
considered them separately in the capacitor equation of 
state. While such a distinction appears desirable from a 
conceptual point of view, implementing it in practice in- 
volves a kind of chicken-and-egg problem. As we have 
shown in Ref. 0, screening and interface bonding are 
strongly interrelated, and it is not obvious how to distin- 
guish cause and effect. At first sight, our approach does 
not appear to provide a solution to the above dilemma, 
as we cast "everything" (chemical bonding and short- 
and long-range Coulomb interactions) into a black box, 
namely, the interface equation of state. However, on 
closer inspection our method does implicitly provide such 
a separation. In fact, by explicitly working as a func- 
tion of a controlled field (here the D field), we auto- 
matically ensure that only the ingredients that are elec- 
trical in nature are included in the interface equations 
of state, Vi{d). The non-polar contributions, which are 
short-ranged and do not have any direct impact on the 
electrostatics, merely enter the definition of the zero of 
the energy, and are therefore implicitly (but rigorously) 
singled out. 

Finally, we would like to comment briefly on our choice 
of D as electrical variable; such a choice is rather con- 
venient, as we have shown in this work. Interestingly, 
using D is frowned upon by some authors in the Landau- 
theory community (Sii^ especially in cases where depolar- 
izing effects are present. A detailed discussion of this is- 
sue would bring us far from the main scope of our work. 
We limit ourselves to noting that, by means of our fixed- 
ly approach, we seek the electronic and structural ground 
state at a given D within a parameter space spanned by 
all the microscopic degrees of freedom (which are implic- 
itly present in our energy functional). This means that 
our description fully accounts for the effects of the "back- 
ground permittivity," of hypothetical competing instabil- 
ities, and of electromechanical couplings. For this reason. 
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it is free from the shortcomings described in Ref. 0. An 
important point to stress is that, in well-behaved cases, 
the electrical equation of state of a given system leads 
to exactly the same description of the physics regard- 
less of which independent variable (electric field £, po- 
larization P, or electric displacement D) is used. This 
point is obvious in linear dielectrics, where P = and 
D = e£. It still holds in non-linear dielectrics that have 
a single energy minimum as a function of the applied 
field £. In more problematic cases, the equation of state 
might become multivalued or even singular, depending 
on the choice of independent variable. In our experience, 
D tends to be a very convenient choice, especially in these 
"difficult" situations. 



C. Relationship to the efFective-Hamiltonian 
approach 

One of the strengths of our approach is the ability to re- 
cast all the complexity of the interface interactions into a 
smooth function of a single electrical variable. This natu- 
rally leads to powerful modeling strategies, as we demon- 
strated in practice in Sections IIIII and IIVI An obvious 
next step would be to combine the interface information 
derived from first principles with higher-level "effective 
Hamiltonian" descriptions of the ferroelectric film,**^'"*^ in 
order to describe phenomena that involve larger length 
and/or time scales (e.g., ferroelectric switching). 

Effective Hamiltonian approaches have been used quite 
intensely in the past few years to investigate size effects in 
ferroelectric nanostructures^^i^i^i^ such as films, wires 
and dots. Generally, the effective Hamiltonian is formu- 
lated and fitted in order to describe bulk behavior, as 
follows. One first identifies a reduced set of local- mode 
variables to describe the amplitudes of the soft ferroelec- 
tric mode and strains within each unit cell. Then a model 
of the energy, written as an expansion in these local- 
mode degrees of freedom, is constructed, and the param- 
eters in the expansion are fitted to a database of first- 
principles calculations. Among the parameters deter- 
mined in this way are some that correspond to short- and 
long-range dipolar interactions. The finite-temperature 
statistical behavior of the system can then be simulated 
using Monte Carlo or molecular-dynamics techniques. 

In order to simulate a nanostructure, the effective 
Hamiltonian description is typically applied without 
modification to a 2D, ID or OD system. This approach 
has allowed for important conceptual advances, for in- 
stance by elucidating the properties of polarization vor- 
tices in nanodisks and nanorodsf^ in such configurations 
the lower dimensionality is largely responsible for the pe- 
culiar behavior. However, in the case of 2D systems such 
as ferroelectric superlattices or thin-film capacitors, an 
important conclusion of our work is that the fine details 
of the interface bonding and electrostatics are crucial to 
determining the overall physical behavior of the system. 
In that sense, the simple abrupt truncation of the dipolar 



interactions^S, which is assumed in the HcS simulations 
discussed above may fall short of faithfully reproduc- 
ing the overall response of a realistic deviceJ- Including 
the interface-specific information in a thin-film effective 
Hamiltonian model would therefore be very desirable. 

For example, for a given interface, one could evalu- 
ate the interface EOS of the HeS at zero temperature 
(T — 0) and compare with the one computed ab-initio 
using the methods described in this work. One could 
then modify the parameters of the Hcs in the vicinity 
of the interface until the T = interface EOS agrees 
with the first-principles one. This would then enable 
one to answer important questions that cannot be di- 
rectly addressed from first principles. For example, what 
is the temperature dependence of the interface equation 
of state? Or, what is the impact of the electrode on the 
stability of the monodomain state versus a polydomain 
one? Finally, by making use of the "locality principle" 
discussed in this work, it would be relatively easy to ana- 
lyze the iJeff results and compare them to the fully first- 
principles values, with significant benefit for both theo- 
ries. To substantiate these arguments, in the following 
we shall briefly discuss two selected Hcs works that are 
particularly relevant in light of our proposed strategy. 

In Ref. [47!, Bin-Omram et al. investigate the impact of 
electrical and mechanical boundary conditions on the po- 
larization and strain of BaTiOs and Pb(Zr,Ti)03 (PZT) 
films. For a BaTiOa film at a 2.0% compressive strain, 
which roughly corresponds to the SrTiOa substrate as- 
sumed in our calculations, the authors of Ref. find a 
spontaneous polarization that increases for thinner films, 
with a value of 0.56 C/ir? at a thickness of 6 unit cells. 
This is qualitatively similar to the effects we discussed in 
Section Hill for Pt/BaTiOs/Pt capacitors, although sig- 
nificantly larger in magnitude. We stress that such an 
enhancement is far from being a systematic property of 
electroded BaTiOa filmsj^. as we suggested above, the 
surface terms in i?eff should be adapted to the specific 
electrode interface on a case-by-case basis. 

In Ref. m the authors report that 2D, ID and OD fer- 
roelectric nanostructures are characterized by "dielectric 
anomalies" in the form of a negative internal suscepti- 
bility x''"*^- It is not unreasonable to think that the 
surface-induced enhancement of the ferroelectric instabil- 
ity, which is built into most H^s thin-film models, may 
be largely responsible for the reported negative in 
the 2D case. Note that x^™*^ is defined there as an av- 
erage over the volume of the film, unlike our definition 
of the local dielectric permittivity e(a;),— which was in- 
troduced in Refs. ^ and'2Q'. This further highlights the 
conceptual advantage of rigorously separating bulk and 
surface/interface effects by performing a local analysis, 
rather than global averages. 

Overall, we believe that the methods developed in this 
work open interesting new avenues for the accurate sim- 
ulation of ferroelectric nanostructures within the H^s 
framework. 
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D. Interpretation of experimental data 

In this section, we ask how one might best make con- 
tact between an experimental set of electrical measure- 
ments on a series of thin-film capacitors of varying thick- 
ness on the one hand, and the analysis tools we developed 
in Sees. IIIII and \TV\ on the other. Of course, it is best if 
the experiment can approach intrinsic conditions insofar 
as possible. For example, the experimental film should 
ideally be in a monodomain state. This may be hard to 
achieve in short circuit, as depolarizing effects might in- 
duce a multidomain state,— but often may be obtained 
by applying a DC bias to the capacitor as in Ref. [S!]. 
Also, the film should be as free as possible from space 
charges arising from charged defects or trapped carriers, 
which may contribute to band-bending effects not con- 
sidered in the theory. 

Our first prediction is that there exists a value of the 
external bias, A{, that yields the same value do of the 
spontaneous electric displacement regardless of the film 
thickness (provided that the interfaces and the film qual- 
ity are similar). Our second prediction is that, around 
this bias value, the inverse capacitance of the films should 
scale linearly with thickness, with the coefficient of pro- 
portionality being directly related to the bulk permittiv- 
ity. Assuming a small range of biases around A( , we can 
discard the third-order coefficients A3 and write 

Viv(d)=A{ + (d-do)Af^ (47) 

with 

Ai^'^^Al+NAl (48) 

We end up with three coefficients that describe the elec- 
trical properties of the capacitors and their dependence 
on thickness. 

To obtain this information, we believe it is best to rep- 
resent the data in a (P, V) plot as in Fig. [TT] rather than 
a {P,£) plot as is usually done. In this way, one ob- 
tains direct visual insight into the existence of a common 
intersection point {Po,Al). Note that the bulk A2 co- 
efficient provides, as a byproduct, useful information on 
the dielectric properties of the film through Eq. as 
it is independent of the specific interface. (In principle, 
it depends only on the applied strain due to epitaxial 
matching, and of course on the operating temperature.) 
By combining Eqs. (|T7|) and we obtain, for the spon- 
taneous polarization of a short-circuited capacitor of a 



given thickness N, 

Po -Po -[y) Al + NA]^- ^^^^ 

Note that, in typical phenomenological models, the in- 
terface is treated as a linear dielectric, which implies 
Al = (io^2- ^2 is related to the interface capacitance 
Ci and to the effective screening length X^g by 

lAi^CY'^^X^s. (50) 

The factor of one-half on the left-hand side relates to 
the assumption of two equivalent interfaces with identical 
electrical properties. 



VI. CONCLUSIONS 

We have developed a comprehensive methodological 
framework for the computation and analysis of ferro- 
electric capacitors with realistic electrodes. Our method 
is based on density-functional theory, and on recently- 
developed techniques for performing calculations at a 
given value of the electric displacement field. By making 
a rigorous separation between the interface and bulk con- 
tributions to the electrical equation of state of a capac- 
itor, we obtain a compact model, of full first-principles 
accuracy, for the electrical (and piezoelectric) response 
as a function of bias potential and thickness. We expect 
these advances to facilitate the comparison of theory with 
experimental data. We also hope that it will stimulate 
a fruitful interaction with other theoretical approaches 
based, e.g., on Landau theory or effective Hamiltonians. 
Application of similar strategies to investigating the in- 
terface coupling between electric polarization, magnetism 
and other structural degrees of freedom (such as octahe- 
dral tilting) that were not considered here are under way. 
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